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There are many software reliability models which try to predict future 
performance of software based on data generated by the debugging pro- 
cess. Unfortunately, the models appear to be unable to account for the ran- 
dom nature of the data [8]. If we debug the same code multiple times and 
use one of the models to make predictions, we observe intolerable vari- 
ance in the resulting reliability predictions. We believe that data replica- 
tion can remove this variance in laboratory type situations and that it is 
less than scientific to talk about validating a software reliability model 
without considering replication. We also believe that data replication may 
prove to be cost effective in the real world, thus our research centers on 
verification of the need for replication and on methodologies for generat- 
ing replicated data in a cost effective manner. 


We have been pursuing the above in the context of the debugging graph 
[81 by simulation and experimentation. Simulation has been done for the 
Basic model and the Log-Poisson model [4] in a manner similar to that 
previously done for the Jelinski-Moranda and Moranda Geometric models 
[9]. That is reasonable values of the parameters have been assigned and 
used to generate simulated data which is then processed by the models in 
order to determine limitations on their accuracy. The experimentation is a 
continuation of that started by White and Harbison. These experiments 
exploit the existing software and program specimens which are in AIR- 
LAB to measure the performance of reliability models. 


Keith Lirette conducted the simulation study of the basic model and has 
written a paper for his Master’s project under our supervision. His con- 
clusions were that the model performed well on all but the smallest sets of 
data, that data replication improved the accuracy of the model and that the 
confidence intervals provided by the model were inaccurate. We have en- 
closed his paper, as well as his data as part of this report but would like 
to add that he is overly generous in his assessment of the models perfor- 



mance without replication. His work on the Log-Poisson model has not yet 
been written up. 

Chris Cowles has conducted an empirical study measuring the performance 
of the Jelinski-Moranda Model and the Moranda Geometric Model using 
data from one of the LIC programs from RTI. He collected replicated data 
along the most probable (approximately) path through the debugging graph 
and fed this data into the two models. His conclusion was that the 
Geometric Model performed well when given accurate replicated data. 
This work used an environment at AIRLAB developed by Richard White 
and is based on code developed by Bernice Becher of NASA. His MS 
project paper and data is' enclosed. 

We plan to continue this investigation of the potential value of replication 
in the context of the debugging graph. The work so far has shown some 
value in this approach but we have not yet completed our work using the 
LIC programs, in that we have been working with a limited part of the de- 
bugging graph for one version of the LIC code. The next stage is to com- 
plete the error graph for this version prior to investigating the other useful 
version of the LIC software available. It is also our hope to apply similar 
techniques to the GCS specimens. 
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ABSTRACT 


Accurate reliability models which predict future failure 
times of software are useful during software development. The 
accuracy of ^models is tested by comparing predictions which use 
the MTTFs of i errors with data from a known i+lst error. 
Evaluations of a model's accuracy are found by tracing a path 
through the error graph of an existing piece of software. It is 
noted that a model's performance may be distorted when an 
abnormally large or small error is used in the sequence of 
errors, and that better results may be obtained by starting the 
model's run after the large dominant errors have been removed. 

It is also noted that to increase the accuracy of predictions, a 
sufficient number of inputs must be used to calculate failure 
rates and MTTFs . 

INTRODUCTION 


In any large software project, it is highly desirable to be 
able to predict the reliability of the software during the 
debugging phase. Only when the software is free (or nearly free 
to within a certain degree of tolerance) of errors is it ready to 
be released. But it is difficult to say exactly when a piece of 
software is nearly free of errors. Latent errors that remain may 
occur later when certain new demands are made on the software 
under untested (but valid) conditions. Reliability predictions 
themselves are at best good estimators of parameters such as the 



number of initial errors or the expected time lapse until the 
next error. 

Mathematical models of software reliability have been 
proposed in the past. [2] [3] [4] These may be used to estimate 
and/or predict software reliability. If a piece of software' is 
reliable enough to allow, say, one error in a million runs to 
pass through to the output, it may be deemed acceptable, ie. , 
good enough for release. The degree of acceptable reliability 
will vary greatly between software packages: for example, one 
error per year in a database application may be acceptable, but 
not for a real-time aviation application. 

The Jelinski / Moranda and the Moranda Geometric De- 
Eutrophication models were used in this project. The JM model 
assumes that the errors found in the software are all the same 
size and have an equal probability of occurring. Further, it 
assumes that there are a finite number (N) of bugs to begin with, 
and that bug fixes are perfect, ie., no new errors are introduced 
when others are corrected. Given a sequence of n time intervals 
between n+i successive errors, an estimate can be obtained for N 
by using maximum likelihood techniques. 

The Moranda Geometric De-Eutrophication model assumes that 
the detection rate will geometrically decrease as each new error 
is found and removed. Estimates can be found for D (the initial 
failure rate) and for k (the rate at which the failure rate 


declines when errors are fixed) . These values can be used to 
calculate an estimate of the failure rate lambda. Like the JM 
model it assumes perfect fixes, uses maximum likelihood 
techniques and uses the sequence of n inter-failure times as 
input. However, unlike the JM model, the Geometric De- 
Eutrophication assumes an infinite number of bugs. 

This project uses an existing piece of software that has 
twelve known errors, corrections available for each error, and an 
oracle constructed from a trusted (gold) version of the software. 
Since errors can be arbitrarily removed from the original 
software, their removal can be done in any desired order. Inputs 
can be run on each of the versions, and each result can be judged 
by the oracle as success or failure. After recording the results 
of testing for each software version, a reliability rate, failure 
rate, and MTTF can be found. If all of the known errors are 
removed one at a time (in some order) , a complete list of MTTFs 
can be obtained for this simulated debugging process. Once 
completed, a sub-sequence of the MTTFs can then be used as input 
to the two models in this project, and predicted MTTFs can be 
obtained. Since the actual MTTF will already have been 
calculated, the models' predicted MTTFs can be compared to the 
real MTTF. Thus, the accuracy of these two models can be 
studied. 

Because the bugs can be removed selectively from this 
candidate software, comparing an estimated MTTF with an actual 



MTTF is possible. These steps are not possible with software in 
the development process (normally no oracle exists at that time) , 
but it is the performance of the software reliability models that 
is being evaluated here. 


DESCRIPTION OF THE SOFTWARE 


The software used in this project is called the LIC (Launch 
Interceptor Conditions) program. It is a Fortran program which 
simulates the testing of certain conditions, and the results of 
these tests are used to decide whether or not to issue a command 
to launch an interceptor missile. The multiple versions of the 
software, the gold version, the control program, are all written 
in VAX FORTRAN, were produced by the Research Triangle Institute 
[1], and exist at the AIRLAB facility at NASA Langley Research 
Center in Hampton, Virginia. The data collection processes were 
run on DEC VAXs using the VMS operating system. 

The existing LIC software was chosen for this project for 
two reasons. First, work on software's error graph (described 
below) had already begun as part of a Master's project by White & 
Harbison[5] in 1990. This project involved work on the top and 
bottom levels of the error graph and further investigation was 
recommended. Second, most of the support apparatus needed for 
the experiment (gold program, control program, random number 
generator, etc.) already existed at NASA Langley. 



Further, there were reasons for originally using the LIC 
program in the White / Harbison project. Not only had the errors 
that existed in the program been identified, but the corrections 
for these bugs had already been written. The versions of the 
software used in the White / Harbison work had the original 
errors left in place but commented out, and the inserted 
corrections were all identified with comment-tags. Thus newly 
required versions of the software used in this project (ie., 
versions with correction combinations not previously required) 
could easily be constructed. 

The naming convention for the different software versions 
used in this project is the same as that used by White / 

Harbison. The individual bugs are identified by numbers 1 
through 12. The number of the errors corrected ("fixed' 1 ) in a 
version are identified in the name of the Fortran program; eg. , 
Fixl235_AND i _6 has five different bugs removed (namely bugs 1, 2, 
3, 5, and 6), and Fixl_T0_9 has the first nine bugs removed. 
Examples of programs using this naming scheme appear in the 
appendices in this paper. The original oracle is based on a 
trusted program called Gold. 

PREVIOUS WORK ON THE LIC PARTIAL ERROR GRAPH 

A complete error graph is a graph with a single top node, a 
single bottom node, and some intermediate levels between them. 

The top node represents the software in its initial state with N 



unknown bugs. The bottom node represents the final version of 
the software with all errors removed. Each horizontal level 
(with multiple nodes) between these single nodes represents a 
stage with some number of bugs corrected: the first level below 
the top node contains all software versions with one error 
corrected, the second level with all combinations of two errors 
corrected, etc. Between the nodes at each level an edge can be 
drawn to each node on the next level that represents a single 
step in the software debugging process. Hence, if there are N 
errors, level one will contain N nodes, each representing a 
single fix. From each node in level one there will be an edge to 
each of the N-l nodes on the second level. Any walk from the top 
node to the bottom node constitutes a path through the error 
graph and represents a possible debugging sequence. 

Each node can be labeled with its reliability measure, which 
is the ratio of successful runs to total tests for the version 
represented by that node. All inputs are randomly generated 
according to the input distribution and a successful run is 
defined as one where the results of that version of the software 
is the same as that of the gold version. A "delta reliability" 
number can be computed that indicates the percentage of change in 
reliability from the previous node. This delta reliability 
indicates the size of a single error within that previous version 
of the software. 


A partial error graph is suitable to display information 


about n known fixes. Since N > n, the partial error graph will 
be smaller and show only paths involving the known fixes. In the 
White / Harbison project, the levels with one bug corrected and 
n-1 bugs corrected were completed. For the one-bug level, the 
initial version of the software was separately fitted with each 
single correction, and the reliability of each such version was 
tested against the gold version. At the (n-l)st level, each 
version was created by installing all but one fix. Again, the 
reliability 'of each of these versions was tested against the gold 
version. 

DESCRIPTION OF THE SOFTWARE TESTING PROCESS 

In order to obtain a sequence of MTTFs with which to judge 
the two models, it was decided to follow a single path through 
the error graph. This path was determined by identifying the 
largest error remaining at each level and correcting it to arrive 
at the node on the next level. This process was repeated until 
all of the known bugs had been removed from the LIC software. 

The determination of the size of the twelve errors was made from 
the first level of the White / Harbison paper. The ordering of 
bugs 1 through 12 by decreasing failure rate is: #1, #2, #3, 

# 5, # 6, # 8, # 7, # 4, # 9, # 10, # 11, # 12. 

It was assumed at each stage of this project that the next 
largest error still to be found along the single path would be 
the same as that in the ordered listing of bugs at level one. 


This decision does not address the question of whether or not the 
removal of bugs one at a time would in any way affect the 
relative size of the remaining bugs. In order to determine the 
largest error of all of those errors left at each stage, it would 
have been necessary at each stage to couple each of the remaining 
errors separately with the corrections done so far along the 
path. Only by trying each possibility could one be sure that the 
identity of the largest error had been found. In this project, 
sixty-six (11 + 10 + ... +1) different software versions would 
have been required. Since each data point required one million 
random inputs and approximately three days of machine time, the 
time required to test along these lines made it prohibitive. It 
is suggested that there may be a need for further study. 

Initially, one million random inputs were done on each 
software version being tested on the VAX and comparisons were 
made against the gold version. The results were categorized by 
the control program as success, abort, or failure; since both 
aborts and failures are disagreements with the results of the 
oracle, they are summed as "Total Failures" in Appendix A. For 
each group of runs, a corresponding second chart in the appendix 
shows the calculated reliability (number of successes divided by 
number of runs) , failure rate (one minus reliability) , and MTTF 
(the inverse of the failure rate) . 

Both reliability models were then checked using the 
empirical data. Appendix B lists the results of the JM model 



using the data from one million runs. Estimates were derived 
from every possible sequence of two or more consecutive 
interfailure times which began with the MTTF for error #1. Also 
used were smaller sequences of consecutive fixes in the amounts 
of 3, 5 and 7 fixes. 

In the JM model the actual value of N must, by definition, 
be larger than n, the number of errors used to predict N, since 
each of the n known errors is one of the N original errors. If, 
for example, the first six inter-time failures are used to find 
N, then the prediction of N includes these six errors. A 
predicted value of N = 20 in this case would indicate that 14 (20 
- 6) errors remain. Also, when using any sub-sequence that does 
not include some of the initial errors, those errors must be 
excluded in determining the number of errors that remain. For 
example, suppose that when using the sub-sequence of MTTFs #5, #6 
and #7 that N is predicted to be ten. This indicates that the 
total error content is ten beginning at error #5, ie., error #5 
is considered to be the first error, since it was the first one 
used in this particular sub-sequence. Previous errors are not 
taken into account here. Since three MTTFs were used to make 
this prediction, seven errors (10 - 3) now remain. 

The results of a typical removal of bugs in the JM model is 
shown in Figure 1 (next page) . The graph approaches the x-axis 
in a step-wise manner; correcting each bug brings the otherwise 
horizontal graph closer by a step toward that axis. It is 


possible for one of two things to happen: If the step-decreases 
were small enough each time, then the graph would not cross the 
axis on the next step and the prediction would be that the 
software still contains some number of errors. Otherwise, if the 
step sizes were large enough so that the next step would force 
the graph to be below the axis, then one would surmise that the 
software was now perfect (ie., free of errors). 


ERROR 

CONTENT 



SEQUENCE OF FAILURE TIMES AS TIME 
INCREASES IN JELINSKI / MORANDA 


FIGURE 1. 


This last case is a common occurrence with the JM 
predictions of N (the number of initial errors in the software) 
in this project. In Appendices B and D listing the results of 
the JM model, the software generally appears to be perfect. 
(Appendix D is the results without the large error #1. More is 
said below about the influence of the size of this error.) In 
Appendix B, because of the large size and subsequent impact of 



correcting error #1, the next correction would be expected to 
produce an improvement on the same scale. It then appears as if 
the software has improved so much as to now be free of errors. 

The data from the one million random inputs was next used on 
the Moranda Geometric De-Eutrophication model. It was decided to 
also run random inputs in the amounts of one, ten, one hundred, 
one thousand, ten thousand and one hundred thousand in order to 
be able to check the robustness of this model on data with 
varying quality. As a hypotheses, it was thought that more and 
more accurate predictions would be produced as the number of runs 
increased. All of these results are listed in Appendix C. As in 
the JM model, predictions are also made without using the 
influential error #1. These are listed in Appendix E. Finally, 
also included is Appendix F, listing the two Pascal programs used 
to run the reliability models. 


RESULTS OF THE MODELS 


The candidate software used in this project contained an 
extremely large and influential first error. The last four or 
five errors could all be classified as very small, with little 
effect on the predictability of the two models. In fact, even at 
runs of 100,000, the last four errors do not appear. Aside from 
these errors, the other corrections seem to improve the software 
in a steady, uniform manner. 


Because of the influence of the first error, computations by 
each of the two models were repeated by first correcting error #1 
and then using the resulting software version as the initial 
version. Hence data and consequently the models' predictions 
were calculated as if this large bug did not exist. In the JM 
model, checked again with the data from one million inputs, only 
the first X consecutive MTTF numbers (from 2 to 11 in this case, 
not 2 to 12 when error #1 was included) were used. Attempting 
the sequences of 3, 5, and 7 again would only shift the answers 
already obtained and listed in Appendix B. (For example, the JM 
result with error #1 on error sequence 4, 5, and 6 would be the 
same as the JM result without error #1 on error sequence 3, 4, 
and 5.) These new results are listed in Appendix D. The data 
from excluding error #1 and using inputs of one, ten, one 
hundred, etc. runs was also checked in the Geometric De- 
Eutrophication model. These new results are listed in Appendix 
E. The discussion of the results of all of the runs follows. 

THE JELINSKI / MORANDA MODEL 

The JM model predictions are not always accurate and often 
indicate perfect software when more errors were known to exist. 
The predictions sometimes vacillate with the correction of a 
single additional error. For example, when using a sequence of 3 
MTTFs (Appendix B) , the predicted N at one point shows 4, then 
perfect twice in a row, then 4, then 18. The number of remaining 
errors in these predictions is, respectively, 1, 0, 0, 1, 15. 



The predicted initial error content of 4, when using corrections 
#5 - #7 or #8 - #10, seems to be good, as this indicates one 
error left. (The LIC program at these points has, respectively, 
five and then two more known errors.) Generally, the JM model is 
not this accurate. 

The overall weakness of this model is possibly the result of 
the existence of the large and small error (s) mentioned above, 
and the relatively small number of known errors (twelve) in the 
software. It should be noted that, in Appendices B and D, when 
sequences of the first X calculated MTTFs are used, the software 
is generally deemed to be free of any more errors. In real 
testing, these are the sequences of inter-time failures most 
likely to be used. 

THE GEOMETRIC DE-EUTROPHICATION MODEL 

The Moranda Geometric De-Eutrophication model performed 
better. Also, as might be expected, increasing the number of 
test cases used to establish the MTTF sequence greatly effects 
the results and the accuracy of the predictions. 

In Appendix C, which includes error #1, the effect of the 
large first error is paramount. Initially for one million runs 
(with only two MTTFs used as input) a MTTF is predicted that is 
5.29 times as great as the actual number derived emperically. 

This ratio of estimated MTTF and actual MTTF declines as more 



data is used; this can be partly accounted for by the uniformity 
of the rate at which the software is improved. The closest this 
model comes to matching the estimated MTTF to the actual MTTF is 
when eight data points are entered (ie., when nine errors have 
been removed) . This ratio is 0.861. Thereafter, near the end of 
the input of available data, the predictions are again skewed 
upward, when data from consecutive very small errors is used in 
the predictions. 

Upon removing the effects of error #1 entirely from the 
initial software, the Geometric De-Eutrophication predictions of 
MTTF (found in Appendix E) improve greatly. Even runs consisting 
of only ten thousand inputs show rather remarkable estimates. 

The more accurate results obtained at one million inputs are used 
to calculate the numbers below. Prior to inserting the 
corrections of the last four small errors, the prediction is 
never wrong by more than 27 percent. (It predicts 127 percent of 
the actual MTTF at 5 data points.) The one poor prediction in 
this range is the .0151 ratio (1159 / 76,923) at 6 inputs. The 
next paragraph explains why it is not accurate in this case. The 
most accurate prediction is actually made with just the first two 
MTTFs : this ratio is just 1.0247 (63.95 / 62.41). 

The De-Eutrophication model at one point seems to "correct 
itself”. Note first, using the failure numbers from Appendix A 
under one million runs, that from Fixl_to_2 to Fixl_to_3 the 
correction is 47 percent ((56,927 - 29,833) / 56,927). The 


succeeding corrections in percents are 46, 60, 52, 37 and then 
99. This 99 percent correction is measured from 1966 failures at 
Fixl23567_and_8 to 13 failures at Fixl_to_8 (ie., (1966 - 13) / 

1966) . The model misses this great correction (at 6 MTTFs 
entered) when it predicts a MTTF of 1,159 and the actual value is 
76,923. However, its next prediction at 7 MTTFs, which now 
includes the seventh value (ie., the failure of only 13 times at 
Fixl_to_8) , is much better. The model has gotten back on track 
with the additional input of the first of the very small errors, 
and the ratio of estimated MTTF and actual MTTF is .8369 (104,612 
/ 125,000). The remaining errors are small and the predictions 
for the rest of the figures are skewed upward a great deal, 
making their ratios very much larger than one. 

The overall stability of the measure of failure rates 
against the gold version improves as the number of randomly 
generated inputs increases. This is not surprising; however, 
even between 100,000 and 1,000,000 runs, the difference in MTTF 
at Fixl2356_And_8 (six errors taken out) is still greater than 
one. (It changes from 319.49 to 322.48.) The other five MTTFs 
before Fixl2356_And_8 are at least stable to within plus or minus 
one. The failure rate of Fixl (which will eventually stabilize 
at about 57.6%) moves from 100% to 40% to 60% from using one, ten 
and then one hundred inputs. Ten thousand runs are required to 
get this to stabilize to within one percentage point. And even 
the last 4 errors only begin to show up at one million runs. 

These changing MTTFs, which are as accurate as the number of runs 



that were made, must have a degrading effect on the prediction of 
the Geometric De-Eutrophication model. 


CONCLUSIONS 


The Geometric De-Eutrophication model seems more adept at 
predicting the time to the next error, if it exists. The 
Jelinski / Moranda model did not perform well. Given the 
software used in this project, with only twelve known errors, the 
predictions were accurate in the Geometric model and might be 
acceptable in a real application. The small number of bugs with 
which to derive data and the results from that data is a credit 
to the Geometric De-Eutrophication model, but may have been a 
factor in hindering the ability of the JM model. 

The elimination of a single abnormally large bug greatly 
improved the Geometric De-Eutrophication predictions. 

Predictions were skewed when this type of error was part of the 
set producing the MTTF data. When this software's large first 
error was removed, and the resulting version was used as the 
initial version, the predicted results were very good. It is 
easy to say that the De-Eutrophication model is more accurate 
when errors of irregular size are ignored, but it is difficult to 
determine when the software is ready to begin being predicted. 

Who is to say when abnormally large errors have ceased and those 
that will follow are decreasing at a uniform rate? If an error 
is classified as large, this does not say anything about the 



errors that follow; they may all be "large 1 ' , but decreasing at a 
uniform rate, as is desired for this model. Perhaps one should 
look for abrupt changes in the rate of increasing MTTFs as 
corrections are made. 

Finally, it is generally true that the greater the number of 
inputs used in order to produce MTTF data, the more accurate the 
predictions will be. This is true throughout the different 
stages of the De-Eutrophication model as shown by the 
increasingly accurate MTTF predictions. The MTTFs used from 
Appendix A are, in some cases, only beginning to stabilize even 
at one million runs. In real-time testing, runs of over one 
million may be desirable. The number of inputs used to predict 
MTTFs might depend on the degree of accuracy needed in the 
application, and on the ability to continue to collect the next 
failure point and introduce it into the data set. At some point, 
comparing past predictions of MTTF (gained from software versions 
with fewer fixes) with those produced using the most current data 
available may help to determine when the predictions are accurate 
enough and the number of inputs used in finding the predictions 
is sufficiently high. 
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Results of random 
runs of 1; 10; 100; 


inputs on all software versions; 

1 , 000 ; 10 , 000 ; 100 , 000 ; 1 , 000 , 000 . 


TOTAL FAILURES 


TOTALS 


Fixl 

Fixl_TO_2 

Fixl_TO_3 

Fixl23_AND_5 

Fixl235_AND_6 

Fixl2356_AND_8 

Fixl23567_AND_8 

Fixl_TO_8 

Fixl_TO_9 

Fixl_TO_10 

Fixl_TO_ll 

Fixl TO 12 


SUCCESS 

0 

0 

0 

0 

0 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


All Gold version results were 100% successful. 


VERSION 

RELIABILITY 

FAILURE 

RATE 

MTTF 

Fixl 

0.0000% 

100.0000% 

1.00 

Fixl And 2 

0.0000% 

100.0000% 

1.00 

Fixl To 3 

0.0000% 

100.0000% 

1.00 

Fixl23 And 5 

0.0000% 

100.0000% 

1.00 

Fixl235 And 6 

0.0000% 

100.0000% 

1.00 

Fixl23 56 And 8 

100.0000% 

0.0000% 

N/A 

Fixl23567 And 8 

100.0000% 

0.0000% 

N/A 

Fixl To 8 

100.0000% 

0.0000% 

N/A 

Fixl To 9 

100.0000% 

0.0000% 

N/A 

Fixl To 10 

100.0000% 

0.0000% 

N/A 

Fixl To 11 

100.0000% 

0.0000% 

N/A 

Fixl To 12 

100.0000% 

0.0000% 

N/A 


RESULTS OF 1 RANDOM INPUT 



SUCCESS 


TOTAL FAILURES 


TOTALS 


Fixl 

6 

4 

10 

Fixl_TO_2 

9 

1 

10 

Fixl_TO_3 

9 

i ? 

• 

10 

Fixl23_AND_5 

9 

1 

10 

Fixl235_AND_6 

9 

1 

10 

Fixl2356_AND_8 

10 

0 

10 

Fixl23567_AND_8 

10 

0 

10 

Fixl_TO_8 

10 

0 

10 

Fixl_TO_9 

10 

0 

10 

Fixl_TO_10 

10 

0 

10 

Fixl_TO_ll 

10 

0 

10 

Fixl_TO_12 

10 

0 

10 

All Gold 

version results 

were 100% successful. 



VERSION 

RELIABILITY 

FAILURE 

RATE 

MTTF 

Fixl 

60.0000% 

40.0000% 

2.50 

Fixl And 2 

90.0000% 

10.0000% 

10.00 

Fixl To 3 

90.0000% 

10.0000% 

10.00 

Fixl23 And 5 

90.0000% 

10.0000% 

10.00 

Fixl235 And 6 

90.0000% 

10.0000% 

10.00 

Fixl2356 And 8 

100.0000% 

0.0000% 

N/A 

Fixl23567 And 8 

100.0000% 

0.0000% 

N/A 

Fixl To 8 

100.0000% 

0.0000% 

N/A 

Fixl To 9 

100.0000% 

0.0000% 

N/A 

Fixl To 10 

100.0000% 

0.0000% 

N/A 

Fixl To 11 

100.0000% 

0.0000% 

N/A 

Fixl_To_12 

100.0000% 

0.0000% 

N/A 


RESULTS OF 10 RANDOM INPUTS 



SUCCESS 

TOTAL FAILURES 

TOTALS 

Fixl 

39 

61 


100 

Fixl_TO_2 

96 

4 


100 

Fixl_TO_3 

97 

3 


100 

Fixl23_AND_5 

99 

1 


100 

Fixl235_AND_6 

99 

1 


100 

Fixl2356_AND_8 

100 

0 


100 

Fixl23567_AND_8 

100 

0 


100 

Fixl_TO_8 

100 

0 


100 

Fixl_TO_9 

100 

0 


100 

Fixl_TO_10 

100 

0 


100 

Fixl_TO_ll 

100 

0 


100 

Fixl_TO_12 

100 

0 


100 

All Gold 

version results were 100% 
FAILURE 

successful 

• 

VERSION 

RELIABILITY 

RATE 

MTTF 


Fixl 

39.0000% 

61.0000% 

1.64 


Fixl And 2 

96.0000% 

4.0000% 

25.00 


Fixl To 3 

97.0000% 

3.0000% 

33.33 


Fixl23 And 5 

99.0000% 

1.0000% 

100.00 


Fixl235 And 6 

99.0000% 

1.0000% 

100.00 


Fixl2356 And 8 

100.0000% 

0.0000% 

N/A 


Fixl23567 And 8 

100.0000% 

0.0000% 

N/A 


Fixl To 8 

100.0000% 

0.0000% 

N/A 


Fixl To 9 

100.0000% 

0.0000% 

N/A 


Fixl TO 10 

100.0000% 

0.0000% 

N/A 


Fixl To 11 

100.0000% 

0.0000% 

N/A 


Fixl_To_12 

100.0000% 

0.0000% 

N/A 



RESULTS OF 100 RANDOM INPUTS 


• 

SUCCESS 

TOTAL FAILURES 

TOTALS 

Fixl 

430 

570 

1,000 

Fixl_T0_2 

947 

53 

1,000 

Fixl_TO_3 

968 

32 

1,000 

Fixl23_AND_5 

989 

11 

1,000 

Fixl235_AND_6 

996 

4 

1,000 

Fixl2356_AND_8 

999 

1 

1,000 

Fixl23567_AND_8 

1,000 

0 

1,000 

Fixl_TO_8 

1,000 

0 

1,000 

Fixl_T0_9 

1,000 

0 

1,000 

Fixl_TO_10 

1,000 

0 

1,000 

Fixl_T0_ll 

1,000 

0 

1,000 

Fixl TO 12 

1,000 

0 

1,000 


All Gold version results were 100% successful. 


VERSION 

RELIABILITY 

FAILURE 

RATE 

MTTF 

Fixl 

43.0000% 

57.0000% 

1.75 

Fixl And 2 

94.7000% 

5.3000% 

18.87 

Fixl To 3 

96.8000% 

3.2000% 

31.25 

Fixl23 And 5 

98.9000% 

1.1000% 

90.91 

Fixl235 And 6 

99.6000% 

0.4000% 

250.00 

Fixl2356 And 8 

99.9000% 

0.1000% 

1,000.00 

Fixl23567 And 8 

100.0000% 

0.0000% 

N/A 

Fixl To 8 

100.0000% 

0.0000% 

N/A 

Fixl To 9 

100.0000% 

0.0000% 

N/A 

Fixl To 10 

100.0000% 

0.0000% 

N/A 

Fixl To 11 

100.0000% 

0.0000% 

N/A 

Fixl To 12 

100.0000% 

0.0000% 

N/A 


RESULTS OF 1,000 RANDOM INPUTS 


o 



SUCCESS 

TOTAL FAILURES 

TOTALS 

• 

Fixl 

4,269 

5,731 

10,000 

Fixl_T0_2 

9,450 

550 

10,000 

Fixl_T0_3 

9,711 

189 

10,000 

Fixl23_AND_5 

9,850 

150 

10,000 

Fixl235_AND_6 

9,945 

55 

10,000 

Fixl2356_AND_8 

9,978 

22 

10,000 

Fixl23567_AND_8 

9,986 

14 

10,000 

Fixl_T0_8 

10,000 

0 

10,000 

Fixl_T0_9 

10,000 

0 

10,000 

Fixl_TO_10 

10,000 

0 

10,000 

Fixl_TO_ll 

10,000 

0 

10,000 

Fixl_TO_12 

10,000 

0 

10,000 

All Gold 

version results 

were 100% successful. 


VERSION 

RELIABILITY 

FAILURE 

RATE 

MTTF 

Fixl 

42.6900% 

57.3100% 

1.74 

Fixl And 2 

94.5000% 

5.5000% 

18.18 

Fixl To 3 

97.1100% 

2.8900% 

34.60 

Fixl23 And 5 

98.5000% 

1.5000% 

66.67 

Fixl235 And 6 

99.4500% 

0.5500% 

181.82 

Fixl2356 And 8 

99.7800% 

0.2200% 

454.55 

Fixl23567 And 8 

99.8600% 

0.1400% 

714.29 

Fixl To 8 

100.0000% 

0.0000% 

N/A 

Fixl To 9 

100.0000% 

0.0000% 

N/A 

Fixl To 10 

100.0000% 

0.0000% 

N/A 

Fixl To 11 

100.0000% 

0.0000% 

N/A 

Fixl_To_12 

100.0000% 

0.0000% 

N/A 


RESULTS OF 10, 

000 RANDOM 

INPUTS 


TOTALS 


' 

SUCCESS 

Fixl 

42,342 

Fixl_TO 2 

94,258 

Fixl_TO_3 

97,011 

Fixl23_AND_5 

98,414 

Fixl235_AND_6 

99,352 

Fixl2356_AND 8 

99,687 

Fixl23567_AND 8 

99,809 

Fixl_TO_8 

99,999 

Fixl_TO_9 

100,000 

Fixl_TO_10 

100,000 

Fixl_T0_ll 

100,000 

Fixl_TO_12 

100,000 


TOTAL FAILURES 


57,658 

100,000 

5,742 

100,000 

2,989 

100,000 

1,586 

100,000 

648 

100,000 

313 

100,000 

191 

100,000 

1 

100,000 

0 

100,000 

0 

100,000 

0 

100,000 

0 

100,000 


All Gold version results were 100% successful. 


VERSION 

RELIABILITY 

FAILURE 

RATE 

MTTF 

Fixl 

42.3420% 

57.6580% 

1.73 

Fixl And 2 

94.2580% 

5.7420% 

17.42 

Fixl TO 3 

97.0110% 

2.9890% 

33.46 

Fixl23 And 5 

98.4140% 

1.5860% 

63.05 

Fixl235 And 6 

99.3520% 

0.6480% 

154.32 

Fixl2356 And. 8 

99.6870% 

0.3130% 

319.49 

Fixl23567 And 8 

99.8090% 

0.1910% 

523.56 

Fixl To 8 

99.9990% 

0.0010% 

100,000.00 

Fixl To 9 

100.0000% 

0.0000% 

N/A 

Fixl To 10 

100.0000% 

0.0000% 

N/A 

Fixl_To_ll 

100.0000% 

0.0000% 

N/A 

Fixl_To 12 

100.0000% 

0.0000% 

N/A 


RESULTS OF 100,000 RANDOM INPUTS 



SUCCESS 


TOTAL FAILURES 


TOTALS 


Fixl 

424,327 

575,673 

1,000,000 

Fixl TO_2 

943 , 073 

56,927 

1,000,000 

Fixl_TO_3 

970,167 

29,833 

1,000,000 

Fixl23_AND_5 

983,977 

16,023 

1,000,000 

Fixl235_AND_6 

993,522 

6,478 

1,000,000 

Fixl2356_AND_8 

996,899 

3,101 

1,000,000 

Fixl23567_AND_8 

998,035 

1,966 

1,000,001 

Fixl_TO_8 

999,988 

13 

1,000,001 

Fixl_TO_9 

999,992 

8 

1,000,000 

Fixl_TO_10 

999,994 

6 

1,000,000 

Fixl_TO_ll 

999,993 

7 

1,000,000 

Fixl_TO_12 

999,994 

6 

1,000,000 


All Gold version results were 100% successful. 


(Runs 

of one million 

one are 

not misprints 



FAILURE 


VERSION 

RELIABILITY 

RATE 

MTTF 

Fixl 

42.4327% 

57.5673% 

1.74, 

Fixl And 2 

94.3073% 

5.6927% 

17.57 

Fixl To 3 

97.0167% 

2.9833% 

33.52 

Fixl23 And 5 

98.3977% 

1.6023% 

62.41 

Fixl235 And. 6 

99.3522% 

0 . 6476% 

154.37 

Fixl2356 And 8 

99.6899% 

0.3101% 

322.48 

Fixl23567 And 8 

99.8034% 

0.1966% 

508.65 

Fixl To 8 

99.9987% 

0.0013% 

76,923.15 

Fixl To 9 

99.9992% 

0.0008% 

125,000.00 

Fixl To 10 

99.9994% 

0.0006% 

166,666.67 

Fixl To 11 

99.9993% 

0.0007% 

142,857.14 

Fixl_To_12 

99.9994% 

0.0006% 

166,666.67 


RESULTS OF 1,000,000 RANDOM INPUTS 


********************* 
* * 

* APPENDIX B * 

* * 
********************* 


Jelinski / Moranda results. 
One million inputs. 
Includes error #1. 



* A * 




- Inter-time 

Predicted initial 


failures entered 

error content 


first 2 

Software now perfect. 


first 3 

Software now perfect. 


first 4 

Software now perfect. 


first 5 

Software now perfect. 


first 6 

Software now perfect. 


first 7 

Software now perfect. 


first 8 

Software now perfect. 


first 9 

Software now perfect. 


first 10 

Software now perfect. 


first 11 

Software now perfect. 


first 12 

25 


Values tested: first X calculated MTTFs. 


Inter-time 

Predicted initial 


failures entered 

error content 


1-3 

Software now perfect. 


2-4 

Software now perfect. 


3-5 

Software now perfect. 


4-6 

Software now perfect. 


5-7 

4 


6-8 

Software now perfect. 


7-9 

Software now perfect. 


8-10 

4 


9' - 11 

18 


10 - 12 

Software now perfect. 


Values tested: 

sequences of 3 MTTF. 


Inter-time 

Predicted initial 


failures entered 

error content 


1-5 

Software now perfect. 


2-6 

Software now perfect. 


3-7 

Software now perfect. 


4-8 

Software now perfect. 


5-9 

Software now perfect. 


6-10 

Software now perfect. 


7-11 

6 


8-12 

10 


Values tested: 

sequences of 5 MTTF. 



********************* 
* * 

* APPENDIX C * 

* * 
********************* 


Moranda Geometric De-Eutrophication results 
Includes error #1. 



n 


Inter-time 
failures entered 

1- 7 

2 - 8 

3- 9 

4- 10 

5- 11 


Values tested: 


Predicted initial 
error content 

Software now perfect. 
Software now perfect. 
Software now perfect. 
Software now perfect. 
Software now perfect. 
8 

sequences of 7 MTTF. 



Number of random inputs: 1 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est . ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.9994 

1.000031 

.999909 

1 

1 

3 

.9994 

1.000061 

.99988 

1 

1 

4 

.9994 

1.000092 

.99985 

1 

1 

5 

.9994 

1.000122 

.99982 

1 

* 


( * Did not fail in next 1 case.) 


Number of random inputs: 10 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

. 25001 

.399989 

.02500 

40 

10 

3 

.49999 

.300007 

.03750 

26.67 

10 

4 

. 68527 

. 231659 

. 05109 

19.57 

10 

5 

.79366 

.194533 

.06126 

16.32 

* 


( * Did not fail in next 10 cases.) 


Number of random inputs: 100 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

. 06558 

.609830 

.00262 

381.24 

33.33 

3 

.22182 

.339927 

.00371 

269.54 

100 

4 

.31629 

.249286 

.00249 

400.82 

100 

5 

.42163 

.173845 

.00232 

431.71 

* 


( * Did not fail in next 100 cases.) 



Number of random inputs: 1,000 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est . ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.09268 

.571610 

.00491 

203.67 

31.25 

3 

.23659 

.376715 

.00499 

200.45 

90.91 

4 

.30970 

.300937 

.00277 

361.21 

250 

5 

.33619 

.273620 

.00118 

851.03 

1,000 

6 

.32740 

.284330 

.00035 

2,855.71 

* 


( * Did not fail in next 1,000 cases.) 


Number of random inputs: 10,000 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.09573 

.574644 

.00527 

189.88 

34.60 

3 

.22426 

.396975 

.00448 

223.35 

66.67 

4 

.33375 

.282887 

.00351 

284.92 

181.82 

5 

.37012 

.249989 

.00174 

575.93 

454.55 

6 

.38440 

.236102 

.00076 

1,312.81 

714.29 

7 

.41040 

.208959 

.00041 

2,440.62 

* 


( * Did not fail in next 10,000 cases.) 


Number of random inputs: 100,000 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.09927 

.578148 

.00570 

175.51 

33.46 

3 

.22743 

.404169 

.00475 

210.32 

63.05 

4 

.33729 

.289201 

.00374 

267.18 

154.32 

5 

.38232 

.248733 

.00203 

492.16 

319.49 

6 

.40918 

.224134 

.00105 

950.65 

523.56 

7 

.43725 

.197623 

.00060 

1,655.95 

100,000.00 

8 

.24635 

.516901 

.000007 

142,606.96 

* 


( * Did not fail in next 100,000 cases.) 


Number of random inputs: 1,000,000 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.09903 

.574725 

.00564 

177.43 

33.52 

3 

.22780 

.400969 

.00474 

210.98 

62.41 

4 

.33948 

.285518 

.00379 

263.69 

154.37 

5 

.38379 

.246333 

.00205 

487.54 

322.48 

6 

.40954 

.223023 

.00105 

950.27 

508.65 

7 

.43872 

.195770 

.00061 

1, 632.95 

76,923.15 

8 

.25709 

.486989 

.0000093 

107,585.69 

125,000.00 

9 

.25502 

.494283 

.0000023 

443,487.55 

166,666.67 

10 

.25783 

.527950 

.0000007 1,459,303.34 

142,857.14 

11 

. 27052 

. 457094 

.0000003 3,853,062.31 

166,666.67 

12 

.28529 

.375454 

.0000001 9,162,777.37 

* 


( * Did not fail in next 1,000,000 cases.) 




********************* 
* * 

* APPENDIX D * 

* * 
********************* 


Jelinski / Moranda results 
One million inputs. 
Excludes error #1. 



- Inter-time 
failures entered 


first 2 
first 3 
first 4 
first 5 
first 6 
first 7 
first 8 
first 9 
first 10 
first 11 


Predicted initial 
error content 


3 

Software now perfect. 
Software now perfect. 
Software now perfect. 
Software now perfect. 
Software now perfect. 
Software now perfect. 
Software now perfect. 
Software now perfect. 
Software now perfect. 


Values tested: first X calculated MTTFs 



********************* 
* * 

* APPENDIX E * 

* * 
********************* 


Moranda Geometric De-Eutrophication results. 
Excludes error #1. 



Number of random inputs: 1 


First n 
calculated 
MTTFs 

k 

• 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.9994 

1.000031 

.99991 

1 

1 

3 

.9994 

1.000061 

.99988 

1 

1 

4 

.9994 

1.000092 

.99985 

1 

* 


( * Did not fail in next 1 case.) 


Number of random inputs: 10 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.99994 

. 100003 

.09999 

10 

10 

3 

.99994 

.100006 

.09999 

10 

10 

4 

.99994 

.100009 

.09998 

10 

* 


( * Did not fail in next 10 cases.) 


Number of random inputs: 100 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.75009 

.040000 

.02250 

44.43 

100 

3 

.49999 

.045002 

.00562 

177.78 

100 

4 

.58555 

.040453 

.00476 

210.28 

* 


( * Did not fail in next 100 cases.) 



Number of random inputs: 1,000 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.60386 

.052993 

.01932 

51.75 

90.91 

3 

.45556 

.057722 

.00546 

183.24 

250.00 

4 

.41577 

.061026 

.00182 

548.37 

1,000.00 

5 

..36817 

.067984 

.00046 

2,174.56 

* 


( * Did not fail in next 1,000 cases.) 


Number of random inputs: 10,000 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.52538 

.055008 

.01518 

65.86 

66.67 

3 

.52220 

.055118 

.00785 

127.41 

181.82 

4 

.46862 

.059044 

.00285 

351.19 

454.55 

5 

.44518 

.061963 

.00108 

922.92 

714.29 

6 

.46032 

.059323 

.00056 

1,771.8 

* 


( * Did not fail in next 10,000 cases.) 


Number of random inputs: 100,000 


First n 
calculated 
MTTFs 

k 

D 

lambda 
(est. ) 

MTTF 
(est. ) 

MTTF 

(actual) 

2 

.52062 

.057406 

.01556 

64.27 

63.05 

3 

.52562 

.057223 

.00831 

120.34 

154.32 

4 

.48742 

.060080 

.00339 

294.90 

319.49 

5 

.47973 

. 061020 

.00155 

645.00 

523.56 

6 

.49352 

.058754 

. 00085 

1,177.97 

100,000.00 

7 

.24843 

.122359 

.00001 

139,949.23 

* 


( * Did not fail in next 100,000 cases.) 


Number of random inputs: 1,000,000 


First n 


calculated 



lambda 

MTTF 

MTTF 

MTTFs 

k 

D 

(est. ) 

(est. ) 

(actual) 

2 

.52416 

.056916 

.01564 

63.95 

62.41 

3 

.53063 

.056680 

.00847 

118.09 

154.37 

4 

.48913 

.059736 

.00342 

292.47 

322.48 

5 

.47973 

.060875 

.00155 

646.53 

508.65 

6 

.49547 

.058309 

.00086 

1,159.16 

76,923.15 

7 

.26027 

.118164 

.0000096 

104,612.75 

125,000.00 

8 

.25758 

.119670 

.0000023 

431,217.65 

166,666.67 

9 

.26661 

. 112449 

.0000008 1,306,476.59 

142,857,14 

10 

.28004 

.100866 

.0000003 3,342,340.95 

166,666.67 

11 

.29774 

.082426 

.0000001 7,442,966.15 

* 


( * Did not fail in next 1,000,000 cases.) 


********************* 
* * 

* APPENDIX F * 

* * 
********************* 


Two pascal programs. 

Jelinski / Moranda reliability model. 
Moranda De-Eutrophication reliability model. 



{ * * -k * *r ** ********************************* 

Christopher Cowles 
ODU Master's Project 
Summer, 1991 

Reliability model by Jelinski & Moranda 

*****************************************} 
program jm (input, output) ; 

const max_try = 1000; {try this many points beyond N} 

data_points = 11; {max points for this model} 

var N, {initial error content} 

m : {number of inter-failure times entered} 

integer ; 

Time : array [1. .data_points] of real; {time intervals between 

successive errors} 

T, {sum of the elements of array Time} 

RHS , {current value of right hand side} 

LHS, {current value of left hand side} 

sigma : {sum from i=l to m of ( i— 1 ) * Time sub_i} 

real; 

lhs_gr eater, {set true if LHS is greater} 

rhs_greater/ {set true if RHS is greater} 

crossed : {initialized to false; ends loop if set to true} 

boolean; 

{****************************************************** 

Function find_sigma 

- This procedure calculates the value of the constant (i-1) 
times Time sub_i. 

*******************************************************}. 
function find_sigma : real; 

var i : integer; 
temp : real; 

begin 

temp : = 0.0; 
for i := 1 to m do 

temp := temp + ((i-l) * (Time[i])); 
find_sigma := temp; 
end; 



{ ****************************************************** 

Procedure Init 

- This procedure reads in the values from the keyboard, 
sets the value of T, sigma, and the 3 booleans. 

*******************************************************j 

procedure init; 
var i : integer; 
begin 

writeln ('How many different time intervals to be entered?'); 
readln (m) ; 

writeln ('Enter the', m:3, ' values, pressing <ENTER> each time.') 
for i := 1 to m do 

readln (Time[i]); {read data into Time array} 

N := m; 

lhs_greater := false; 
rhs_greater := false; 
crossed := false; 
sigma := find_sigma; 

T : = 0.0; 

for i := 1 to m do 

T := T + Time[i]; {sum the elements of array Time} 

end; 

{**********************************£******************* 

Function get_left_side 

- This procedure calculates the value of the left hand 
side of the equation. 

*******************************************************} 
function get_left_side : real; 

var i : integer; 
temp : real; 

begin 

temp : = 0.0; 

for i := 1 to m do 

temp := temp + (1/ (N-(i-l) ) ) ; 

get_left_side := temp; 
end; 



^* ***************************************************** 
Function get_right_side 

- This procedure calculates the value of the right hand 
side of the equation. 

*******************************************************y 
function get_right_side : real; 

begin 

get_right_side := m/ (N - ( ( 1/T) * sigma)); 
end; 

{****************************************************** 
Procedure Compare 

- This procedure determines if the two values (RHS & LHS) 
have crossed in value. 

*******************************************************y 
procedure compare; 

begin 

if rhs_gr eater and (LHS >= RHS) 
then crossed := true 

else if lhs_greater and (RHS >= LHS) 
then crossed := true; 

if crossed then writeln ('The values have crossed at ', N:5) ; 
end; 

{************************************** **************** 
Procedure do_f irst_compare 

- This procedure does the first comparison of LHS & RHS 
at the first value of N. 

******************************************************* j. 

procedure do_f irst_compare; 
begin 

LHS := get_left_side; 

RHS := get_right_side; 
if LHS > RHS then lhs_greater := true 
else if RHS > LHS then rhs_greater := true 
else begin 

crossed := true; 

writeln ('Values are initially equal at ' , 
end; 

end; 


N:4) ; 


^********************** 

MAIN MODULE 

***********************} 

begin 

init; 

do f irst_compare; 

while (N < max_try) and (not (crossed) ) do 
begin 

N := N + 1; 

LHS := get_left_side; 

RHS : = get_right_side; 
compare ; 
end; 

if not (crossed) 

then writeln ('Values never cross up to max_try:5); 
end. 


Pascal program for Jelinski / Moranda model 



{**************************************** 
Christopher Cowles 
ODU Master's Project 
Summer, 1991 

Geometric De-Eutrophication reliability 
model by Moranda 

*****************************************j 

program de_eutrophication (input, output) ; 
const 

step = 0.0001; {get this accurate} 

data_points = 11; 

max_right = 1.0; {maximum value of k} 


var 

m : integer; {successive number of failures} 


RHS , 

LHS, 

test_point, 

left, 

right: 

real; 


{current value of right hand side} 
{current value of left hand side} 
{test this valtie, between 0 and 1} 
{value of left border} 

{value of right border} 


LHS_greater_on_lef t , 
RHS_greater_on_r ight , 
LHS_greater_on_r ight , 
RHS_greater_on_lef t : 
boolean; 


Time : array [1. .data_points] of real; 

{********************************£********************* 
Function power 

- This function calculates x to the n power. 

*******************************************************j 

function power (x, n : real) : real; 
begin 

power := exp (n * (ln(x))); 
end; 


^****************************************************** 
Function get test_point 

- This function sends back a point midway between the 
current values of left and right. 
******************************************************* 

function get_test_point : real; 


begin 

get_test_point := ((right - left)/2) + left; 
end; 

{****************************************************** 

Procedure reassign_boundary 

- This procedure determines if the left or right hand side should be 
set equal to the current test_point. Boundary should not cross 
over the point where the LHS and RHS themselves cross between 0 
and 1. 

******************************************************* y 

procedure reass ign_boundary; 

begin 

if LHS_greater_on_left and (LHS > RHS) 
then left := test_point 
else if RHS_greater_on_left and (RHS > LHS) 
then left := test_point 
else if LHS_greater_on_right and (LHS > RHS) 
then right := test_point 
else if RHS_greater_on_right and (RHS > LHS) 
then right := test point 

else writeln ('ERROR Tn reass ign_boundary at ' , test_point:5:4) ; 

{this last stmt, is a flag} 

end; 

{ **************************************** ************** 

Procedure Evaluate 

- This procedure evaluates LHS and RHS with the value 
of k (between 0 and 1) that is sent to it. 

*******************************************************} 
procedure evaluate (k:real) ; 

var i : integer; 
numerator , 
denominator : real; 


begin 

LHS := (m+1) / 2; 
numerator : = o . o ; 
for i := 1 to m do 

numerator := numerator + (i * (power (k,i)) * Time[i]); 

denominator := 0.0; 
for i := 1 to m do 

denominator := denominator + ( (power (k, i) ) * Time[i]); 
RHS := numerator/ denominator ; 
end; 


{****************************************************** 

Procedure Init 

- This procedure reads in the values from the keyboard, 
sets the value of T, sigma, and the 3 booleans. 

******************************************************* j. 

procedure init; 
var i : integer ; 
begin 

write In ('How many different values to be entered?'); 
readln (m) ; 

writeln ('Enter the', m:3, ' values, pressing <RETURN> each time.') 
for i := 1 to m do 
readln (Time [ i ] ) ; 

LHS_greater_on_left := false; 

RHS_greater_on_right := false; 

LHS_greater_on_right := false; 

RHS_greater_on_left := false; 

left := step; {start at min value of k) 

right := max_right; {start at max value of k} 

evaluate (left) ; 

if LHS > RHS then LHS_greater_on_left := true 
else RHS_greater_on_left := true; 
evaluate (right) ; 

if RHS > LHS then RHS_greater_on_right := true 
else LHS_greater_on_right := true; 
end ; 

{****************************************************** 

Procedure do_report 

- This procedure prints out the desired results. 

*******************************************************} 
procedure do_report; 

var D, 

lambda , 
denominator 
: real; 
i : integer; 

begin 

writeln ('The value of k is ', test_point: 6: 5) ; 
denominator := 0.0; 
for i := 1 to m do 

denominator := denominator + ( (power (test_point, (i-1))) * 

Time[i] ) ; 

D := m / denominator; 
writeln ('The value of D is ' 


{failure rate) 

{denominator of fraction to find D} 


/ 


D : 8 : 6 ) ; 


lambda. := (D * power (test_point, m) ) ; 
writeln ('The failure rate lambda is estimated at ' , 
lambda :6:7) ;{6:5) ;} 

writeln ('The MTTF is estimated at ', (1/lambda) : 10: 2) ; 
end ; 

{********************** 

MAIN MODULE 

*************** ********y 

begin 

init; 

while ( (right - left) > step) do 
begin 

test_point := get test_point; {get new mid-point} 
evaluate (test_poxnt) ; {evaluate LHS and RHS} 

reassign_boundary; {set either left or right to 

test_point} 

end; 

do_report ; 
end. 


Pascal program for Moranda Geometric De-Eutrophication model 
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Simulation Studies of The Basic Execution Time Model 


by Keith J. Lirette 


Abstract 

The Basic Execution Time software reliability model was tested 
for reliability parameter estimation accuracy on simulated data 
from the distribution it is based on. The model performed well on 
all but the smallest sets of test data generated. Data replication 
was tested and proved to increase the accuracy of parameter 
estimations. Confidence intervals generated proved to be 
ineffective in providing an accurate range for estimates. 

0. Introduction 

Software reliability is the probability that a computer 
program will operate failure-free for a particular period of time 
in a particular environment [1] . A large portion of software 
reliability literature is aimed at reliability estimation. 
Mathematical software reliability models have been developed which 
attempt to provide estimations of software reliability via 
statistical methods. Generally, these models use failure time data 
of the software in question and make estimations of the total fault 
content of the software and the current software failure rate. The 
difference between models lies in the assumptions made about the 
probability distribution of program failure times and the relative 
effect each fault has on the total failure rate of the software. 

This paper is devoted to the study of one software reliability 
model. The focus of the study is to determine how well the model 
estimates reliability parameters from data generated according to 
the assumptions made by the model (distribution of faults, ...). 
If the model makes accurate estimates, then it is up to the 
potential model user to determine if the assumptions made by the 
model are reasonable for the software in question. If the model 
fails to make accurate predictions, then the internal consistency 
of the model should be questioned. 

The paper continues on to study the effect of data replication 
on estimates and the use of confidence intervals as a way of 
providing a probable range rather than a single point estimate. 

1 . The Model 

The model studied is the Basic Execution Time Model developed 
by John Musa. According to Musa, this model "is the most 
thoroughly developed and has been most widely applied to actual 
projects." The model is characterized by an exponential failure 
intensity (with respect to time) and a non-homogeneous poisson 
failure distribution [1] . 



The model is based on the assumptions [1] : 

1) failures occur as a random process, which is a non- 
homogeneous poisson process, 

2) a fixed failure intensity per fault (i.e. all faults 
contribute equally to the overall failure intensity) , 

3) fault corrections are immediate upon detection and 
perfect . 

The failure intensity as a function of the number of faults 
removed is given by the function [1] : 


1(H) 


= In 


1 " 


JL 

Yo 


H equals failures experienced, 1 0 is the initial 
failure intensity, and y 0 total faults in the system. 


The graph below shows the effect of removing faults on the failure 
intensity [1] . 



The failure intensity as a function of time is given by the 
exponential equation [1] : 

l(t) = l 0 exp[-^t] 
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For generalization purposes, Musa reparameterizes this equation to 
[ 1 ] 

X(t)=p 0 P 1 exp(-P 1 t) 

where p o = y 0 , estimated total failures, and 
Pi = — , estimated failure rate per fault. 


The graph below shows an example of passing execution time on the 
failure intensity [1] . 


Faitire Intensity 


Initial 

Failure 
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Execution Time 


Parameter estimate equations are derived by the use of maximum 
likelihood techniques. Maximum likelihood can be used when the 
underlying distribution of the data is known. Maximum likelihood 
yields parameter estimates that make it most likely that the 
collected data occurred [1] . The mathematics is detailed in 
Appendix F. The maximum likelihood equations for B 0 and Bj. are 
given by [1] : 


0 1 - exp (Pi t a ) 

and [1] 

_ m a^~a _ V"' £ _ q 

Pi exp (Pita) -1 1 

where m e is the number of faults found and fixed, 

t e is the time to end of testing, and 

tj/s are the times to failure for each fault found. 
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2. Data Simulation and Estimation 


The equations described in Section 1 provide the theory to 
base an estimation program upon. The first step in testing the 
model is the simulation of failure data. This was accomplished by 
using the distribution function of faults [1] : 

C+e«uH 

- [ J. (Jrt dx 

Fit) = 1 - e *— 

where t = time, and tsum is the sum of all 
previous failure times . 


The model assumes that failures occur randomly. Thus we can 
randomly generate values between 0 and 1 for F(t) and solve the 
right hand side of the distribution function for t (which is a 
simulated interfailure time) . Repeating this, we can generate a 
set of failure times. The time variable t is solved for as 


-Yo 


In 




In (1 - r) 


Yo e 


■Uo/To> Csum 


where, r is the randomly generated value for F(t) . 

Implementing program solutions for the maximum likelihood 
equations covered in Section 1, estimates for the B 0 and B x 
parameters can be found for generated failure time sets. The 
program "BASIC. P" contained in Appendix G takes inputs for the 
total number of program faults, the failure rate per fault, the 
number of failure times to simulate, and the number of failure time 
sets to generate (50 was used throughout) . The program generates 
the desired number of failure time sets and makes estimates for B 0 
and B t . 

3. Modeling Results 

Parameters were chosen to produce failure data with. Three 
different values for the total number of program faults were 
used : 50, 100, and 200. A series of trials were run with each 
total fault count. The number of failure times generated started 
at 10% of the total fault count and was increased in 10% increments 
until 90% of the total program faults were generated (e.g. for 100 
total program faults, 10, 20, ..., 90 trial sizes were used) . Each 
trial consisted of estimates based on each of 50 failure time set 
generations (50 for 10 failures, 50 for 20 failures, etc.). The 
same per fault failure rate was used throughout, 0.001. Scatter 
plots were produced showing the difference between the estimates 
and the actual parameter value. Additionally, summary graphs were 
generated to allow for the detection of trends. 

Appendix A contains the scatter plots and summary graphs for 
each of the trials. The scatter plots show the percent difference 
of each of the 50 calculated estimates of total program faults, B 0 , 
from the actual parameter value. An entry in a particular column 
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indicates that an estimate was within plus or minus 5% of that 
amount off of the actual parameter (e.g. an entry in the 10 column 
indicates that a particular estimate was between 5% and 15% off 
from the actual parameter) . The summary graphs show the percentage 
off of the mean of the 50 calculated estimates from the actual 
parameter for each of the trials run for a particular total fault 
count . 

For purpose of analysis I needed to define a measure of 
acceptability for the estimates. I chose to use this criteria: a 
set of estimates are acceptable if a significant majority of the 
estimates (2/3) are within plus or minus 30% of the actual 
parameter. The closer the estimates to the actual parameter, the 
better the estimates. This measure was arbitrarily chosen, but is 
illustrative of other possible criteria. 

Using the criteria defined above, only when 90% of the total 
faults were generated was an acceptable set of estimates calculated 
for the 50 total fault size trials. The 100 fault trials reached 
acceptability at 60% of total faults generated, and estimates 
continued to improve with greater number of faults generated. The 
200 fault trials also reached acceptability at 60% of total faults 
generated. 

The 50 total fault size appears to provide too small of an 
event space to allow accurate estimation. Comparison of the space 
sizes for the 100 and 200 fault trials reveals that the estimates 
of the 200 fault case are more accurate. It is hypothesized that 
this is due to the larger number of failure times generated. For 
example, a 50% trial size for 100 total faults is 50 failure times, 
where a 50% trial size for 200 total faults is 100 failure times. 
Note that the expected time to find 50% of the faults is equal for 
all trial sizes (see t e mean in summary plots Appendix A) . With 
the same per fault failure rate, a program with 200 faults will 
fail at double the frequency of a program with 100 faults. Thus, 
twice the number of failures are experienced in the same span of 
time . 


It was also observed that with the same number of failure 
times generated, estimation for the smaller total fault count is 
more accurate. It is hypothesized that this is due to fact that 
removal of faults in the program with the smaller fault count has 
a larger effect on the residual program failure intensity, causing 
greater times between failure, and allowing for better statistical 
analysis . 

Sample tests were performed with different per fault failure 
rates. The results scaled perfectly. That is, if 0.002 was used 
instead of 0.001, estimates for B x were exactly double, and 
estimates for B 0 were exactly the same. 
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4. Data Replication 

Nagel [2] introduced the idea of data replication as a means 
to more accurately study the assumption that all faults have the 
same failure rate. Shen and Wilson [3] used data replication to 
test its effects on the accuracy of parameter estimation by the 
Jelinski-Moranda and Moranda Geometric reliability models. I also 
used it to test its effect on the accuracy of reliability parameter 
estimation. Instead of producing a single set of failure data, r 
replicates are made. This is synonymous with debugging the same 
program from its initial version r times. The corresponding 
failure times are averaged to obtain a single set of failure times. 
This average of the replicates is then used to estimate the 
reliability parameters. 

Trials were run, each consisting of 50 separate failure time 
generations as before, but with each time set consisting of 2, 3, 
4, 5, 10, 20, 30, ..., or 100 replications. Appendix B contains 
scatter plots and summary graphs like those in Appendix A, but for 
the 30 replicate data trials. Comparison of the scatter plots for 
no replication in Appendix A with those with 30 replications in 
Appendix B, displays the sharp improvement in estimation. For the 
50 total fault trials, estimations become "acceptable" with 70% of 
faults generated and 30 replications versus 90% with no 
replication. For the 100 total fault trials, acceptability is 
reached with 50% of faults generated versus 60% without 
replication. And, for the 200 total fault trials, acceptability is 
reached with 40% of faults generated versus 60% without 
replication . 

Appendix C contains summary graphs which show the effect of 
increasing the number of replications for a given number of failure 
times with each total fault count used. Each Appendix C graph 
shows the change in the B 0 parameter estimation for one trial size 
with the different number of data replications used. Study of the 
summary graphs shows that for the 50 total fault trials, parameter 
estimation improvement begins at 25 failure times produced and 20 
data replications. Review of all the summary graphs indicates that 
parameter estimation improvement is minimal after 30 replications . 
The summary graphs behave asymptotically starting at 30 
replications and continuing onto 100 replications. Replication 
does provide greater estimation improvement for smaller total fault 
counts. This is because smaller total fault counts provide less 
accurate estimations and thus have greater room for improvement. 
Still, with the same amount of replication and percentage of total 
faults produced, larger total fault counts yield more accurate 
estimates than smaller total fault counts. 

What does replication provide? It removes some of the 
randomness of fault occurrence and provides a truer picture of the 
distribution of program faults. The figures contained in Appendix 
D, illustrate the gains from replication. It is shown that by 
increasing replication, fewer faults need to be located to achieve 
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the same accuracy of estimation as for fewer replications. 


5. Confidence Intervals 

Confidence intervals allow for generation of a range of 
parameter estimates that explain the data collected with some level 
of confidence. Musa provides the equations required to generate 
confidence intervals for B 0 and B x . Once a confidence interval is 
obtained for B lf the endpoints of its confidence interval can be 
substituted into the equation for B 0 to obtain a corresponding 
confidence interval for B 0 [1] . 

The upper /lower limits of an approximate 
100 (1-a) percent Cl for P x are given by 

■^l-a/2 






where ^. a / 2 is the appropriate normal deviate, and 


X(P L ) is the expected or Fisher information given by: 
1 _ t|exp(P 1 t a ) 
pj [exp (Pi.t a ) -1] : 


KPi) =m e 


One problem encountered was that the offset for the confidence 
interval around B x was often larger than the estimate for B t . Thus, 
the lower bound of the B x confidence interval was truncated at 0. 
But, this value was needed to calculate the upper bound of the 
confidence interval around B 0 . Zero could not be used in the 
equation to solve for B 0 (this would cause a zero divide) . 
Additionally, testing showed that confidence intervals calculated 
with larger percents of the total faults generated (70, 80, 90%) 
did not provide the required confidence. 

Appendix E contains summary graphs of 50, 55 , ..., 95 percent 
confidence interval tests for the 3 total fault count sizes 
(50,100,200). Data was generated in the same manner as for the 
point estimation graphs: starting at 10% of the total faults and 
increasing by 10% until 90% of total faults were generated. 200 
runs were made on each trial size and an appropriate confidence 
interval calculated for each. For each of the 200 confidence 
intervals calculated, it was determined whether the actual 
parameter value used to generate the estimate fell within it. The 
summary graphs show the percent of the 200 runs whose confidence 
interval contained the actual parameter value. 

The accuracy of the confidence intervals, in relation to how 
many contained the actual parameter value, remained fairly constant 
between the different size total fault count trials. With smaller 
numbers of faults generated, parameter estimates were less accurate 
and confidence intervals were extremely large and useless. As 
parameter estimates improved with a greater number of total faults 
generated, confidence intervals tightened to a useable size. But, 
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these confidence intervals did not contain the required percentage 
of the actual parameters, and thus do not provide the expected 
confidence. As seen in the scatter plots of B 0 in Appendix A, the 
"acceptable" estimates tended to be greater than the actual 
parameter. Thus, confidence intervals around these greater values 
tended to miss high (i.e. the lower bound of the confidence 
interval was greater than the actual parameter) . 

6. Conclusion 

The Basic Execution Time Model generally performed well on 
test data generated from the distribution it is based on. The 
exception was when a smaller number (50) of program faults were 
used, estimations were less accurate. Other total fault sizes 
performed well when at least 60% of the 'total faults were 
generated. Replication was found to be a means of increasing the 
accuracy of parameter estimations. With replication fewer faults 
needed to be found to get estimations as accurate as those 
estimated without replication. Finally, confidence intervals 
proved to be ineffective in providing a probable range for the 
estimated parameters. 
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Appendix A 


This appendix contains two types of graphs showing the modeling results for 
single repetitions. The first graphs are scatter plot of trialB for the 
different total fault counts used and number of failures generated. Each scatter 
plot graph shows how much the each of the 50 generated estimates for Bq, total 
program faults, differ (within plus or minus 5%) from the actual parameter 
Gammaj. n indicates the number of failures times produced, and phi is the per 
fault failure rate used. 

The second graphs are summary graphs of the scatter plot graphs. Each 
summary graph shows how much the mean of the 50 calculated estimates for 
differs from the actual parameter Gamma^ for the different number of failures 
generated. 
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This appendix contains two types of graphs showing the modeling results for 
thirty replications. The first graphs are scatter plot of trials for the 
different total fault counts used and number of failures generated with 30 
replications. Each scatter plot graph shows how much each of the 50 generated 
estimates for Bg, total program faults, differ (within plus or minus 5%) from the 
actual parameter Gamma., n indicates the number of failures times produced, r 
indicates the number or replications (30), and phi is the per fault failure rate 
used. 


The second graphs are summary graphs of the scatter plot graphs. Each 
summary graph shows how much the mean of the 50 calculated estimates for 
differs from the actual parameter Gammag for the different number of failures 
generated using 30 repetitions. 
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Summary Graph of % Mean Bq Estimate differs from Gamma^ 
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Appendix C 


Appendi:; C contains summary graphs which show the affect of increasing 
the number of replications for a given number of failure times with each total 
fault count used. An Appendi:; C summary graph shows the change in the mean of 
50 B 0 parameter estimations for one trial size with the different number of 
data replications used. 
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Summary Graph of 7. Mean Bq Esp-oqate differs from Gamma^ 
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' Summary Graph of */. Mean Bq Est^~Nte differs from Gamm^, 

Gamma^ = 100 Phi = 0.001 Source : 100qrp70 

n = 70 Output : 100grp70.gph 
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' Summa’ry Graph o-f 7, Mean Bq Est/^'te differs -from Gamm^, 




Gamma q = 100 F'hi = 0.001 
n = 80 
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Summary Graph o-f ’/. Mean B Q Esf'^ate differs -from Gamm^-, 
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Gamma = 100 Phi = 0.001 
n = 90 


Source : 100grp90 
Output : 100grp90.gph 
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' Summary Graph of ’/. Mean Est/’^te differs from Gamms^ 


Gamma f) = 200 Phi = 0.001 

Source : 200grpl0 
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' Summa'ry Graph o-f ’/. Mean B 0 Est/^te differs -from Gamm^-, 
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Summary Graph o-f 7. Mean Bq Estr^te differs -from Gamm^ 


Gamma q = 200 F'hi 
n = 60 
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Summary Graph of ’/. Mean B., Est'^Nte differs from Gamm^ 




Gamma q = 200 F'hi - 0.001 
n = SO 
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Summary Graph of /. Mean Bq Estimate differs from Gamm^ 

a Q = 200 Phi = 0.001 Source : 200grp50 

n = 100 Output : 200grp50.gph 
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Summary Graph of 7. Mean Estimate differs -from Gammas 

Gamma 0 = 200 F'hi = 0.001 Source : 200grp60 

n = 120 * Output : 200grp60.qph 
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Appendix D 


The figures contained in Appendix D, illustrate the gains from replication. 
For each of the total fault counts used (50, 100, 200) , the number of replicates 
required for the B 0 estimate to be specific %s off from the actual parameter with 
the different number of failure times generated are charted. 


Example. On page D-l, the entry '35' in the 50% column, 4 row, indicates that 
for 35 failure times simulated and 4 replicates used, the mean of the 50 
estimates calculated is 50% off from the actual parameter value. 



50 Total Program Faults: Grid entries are Number of Failures needed to be 

generated for the number of Replicates indicated to be %off indicated. 
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100 Total Program Faults: Grid entries are Number of Failures needed to be 

generated for the number of Replicates indicated to be %off indicated. 

% Mean of 50 B 0 Estimates Differs from Actual Parameter 


replicates 

50% 

40% 

30% 

20% 

10% 

5% 

4% 

3% 

2% 

1% 

0% 

1 




70 

80/90 







2 



70 









3 



60 


70 


90 





4 











% 

5 



50 

60 








10 





60 

80 





90 

20 


40 


50 








30 







80 





40 

40 








-90 



50 












60 








-90 




70 






-90 


80 




80 












90 












100 













* indicates that estimate is less than actual parameter, otherwise 

greater. 


D-2 






















'"'N 


200 Total Program Faults: Grid entries are Number of Failures needed to be generated for 

the number of Replicates indicated to be %off indicated. 

% Mean of 50 B 0 Estimates Differs from Actual Parameter 
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This Appendix contains graphs showing the results of Confidence Interval 
testing. 50, 55, 60, — , 95% Confidence Intervals were calculated for 200 
parameter estimates for each total fault count with 10, 20, ...90% of total 
program faults generated. Each graph column shows the percent of the 200 
calculated confidence intervals which contained the actual parameter value. For 
example on page El-1, the in row 85 column 20 indicates that 85% of the the 
confidence intervals calculated using 20 simulated failure times contained the 
actual Bj parameter value used to simulate the failure data. And, the @ in row 
65 column 20 indicates that 65% of the confidence intervals calculated using 20 
simulated failure times contained the actual Bj parameter value used to simulate 
the failure data. 
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Confidence Interval Equations 
Failure Time Generation Equations 



n ^ 

EQUATIONS DEFINING BASIC MODEL 


1. The Failure Intensity Equation of the model is given by: 

X(t) = X 0 exp[-— t] 

Yo 

where t = time, X 0 = estimated initial failure intensity 
and Yo = estimated total number of software faults. 


2 . For generality, the Failure Intensity Equation is reparameterized as 
X{t) =P 0 P 1 exp(-p i t) 

where P Q = Yo» estimated total failures, and 

p x = — , estimated failure rate per fault. 

Yo 

3. The Likelihood Equation is given by: 

L(P) = P[T a ^ > t g \ T„' = t m ] 

where m g equals the no. of faults found/ fixed. 
t m equals time of last fault. 
t ± equals total time of ith failure. 
t e equals time at end of testing. 

T t is the random variable for time of ith failure. 
f(t i \t 1 _ 1 ) is the conitional density function of (T lf . . .,T m ) 

P [ T m ^ > t g | T mt = t a# ] is the conditional distribution of T a ^ 

The Conditional Density Function is given by : 
fftjti-i) = l(t 1 )exp{-[ji(t 1 ) - H ( ) ] } 


where the Mean Value Function (mean faults by time t) is given by: 
H (t) = p 0 [1 - expt-Pit) ] 


The Conditional Distribution is given by: 

\T a =t„] = exp {-[|i(t.) - 

4. Substitution into the Likelihood Equation yields: 
L(p> = [rh(ti)] exp [~|i c t„) ] 

5. The logarithm of the Likelihood Equation is: 

m. 

In L(p) =£)lnl(t i ) - jiCfc,) 

l-i 


F-l 



6. Thus, the Maximum Likelihood Equation becomes: 


£ i axit ± ) 

k Mti) ap* 


a \i ( t 8 ) 

ap* 


= 0, k = 0, 




7. It can be verified by direct substitution into the previous equation with 
k = 0 that all scaled Mean Value Functions yield the result: 

which can be rewritten as 
ft 

0 1 - exp (Pi.t e ) 

% 

8 . 

The probability observing m a failures by time t a is given by 
Tuft ) 1 "• 

P[M{t e ) - m g ] = i exp[-ii(t„)] 

m g \ 


9. After appropriate substitutions, the Conditional Likelihood Equation becomes 


I W e) = 


1-1 




10 . 

The Conditional Likelihood of p t given m 0 failures is: 


L(Pi | m 9 ) 


/n a , Pi*exp(-p 1 V t ± ) 
1=1 

[1 - expC-p^)]®* 


11 . 

Using the In of the maximum conditional likelihood eq. , the 
point estimate for P x is given as: 

£ - - T Ci = o 

Pi exptPjt,,)-! hi 
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Confidence Interval Equations 


12. 


The MLE of P k is asymptotically normally distributed with mean 


variance 


I(Pk> 

information given by: 


where I(p k ) is the expected, Fisher 


J(PJ = E 


yinup^yp) 

ap* 


13. Thus we can write, 
N( 0,1) 


P*-P* 


N 


14. 

The upper/lower limits of an approximate 100 (1-a) percent 
Confidence Interval for p k are given by 

P*± 

where K 1 . a/2 is the appropriate normal deviate. 


15. 


U Pi) 



t«exp(P 1 t g ) 
[exp(P 1 t a ) -l] 2 


Substituting the confidence interval values for B x into Equation 7, a 
corresponding confidence interval is obtained for B„ . 
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Time Generation Equationb 

16. The Distribution of Faults as a function of time is given as 

c*c*ub 

- f 1 (jc) dx 

F(t) - 1 - e ““ 


The simplification of the integral is found to be: 

£*tsum C+esum 


taiun 


f X(x)dx= f X 0 e Uo/To>JC dx 

•n tsum 

- [-Y 0 e' (VT,,x + c]|^T 

_ J_Y Q g-<Vro> te» taunt) + c _ (_Y o g-(VTo)e*un + 


= Yo e 


-<A 0 /Yo> ts 


“tl 


(^o/YoJ * 


17. Substituting the integrated exponent we get: 
F(t) = 1 - e T ° a (1 a 1 


18. Randomly generate values for F(t), r, and solve for t. 
1 - r = e~ 




19. Take the natural logarithm of both sides and simplify. 
ln(l - r) = -Y 0 e- (V,olt! “" d-e' 11 "'} 


20. Rearrange algebraically. 
-ln(l - r) _ 

v --(Vro)tsum 

Yo® 


21. Rearrange algebraically. 

1 +. ln(l - r) = e -a 0 /Y 0 )c 

„ -(X 0 / 1o ) csum 

To® 


22. Take the natural logarithm of both sides. 

t 




23. Rearrange algebraically. 


3l In (l * 


l 1 


Yo e 


“C^o/Tq) Caum 


) 
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Appendix G 


— CodeP 

This Appendix contains the source code of programs used to calculate 
ailure times, parameter estimates, confidence intervals, and produce graphs 
sed. 


APPENDIX G INDEX 


G-l Basic. p : Failure Times and Estimation 

G-6 Confidence. p : Confidence Intervals 

G-17 Graph.p : Summary Graphs 

G-22 Histogram. p : Scatter Plots 



BAS r V P 


This program generates -failure data using the distribution o-f the Basic Execution Time Model. 
It then calculates estimates -for the reliability parameters Bq and Bj and writes these to an 
output -file. The program reads input from a text file. Inputs are the total program faults 
f gamma_0) , failure intensity per fault (phi), the total number of faults to generate (Me), the 
'umber of runs to perform for each failure set size (runs) , and the number of data repititions 
(repeats). The output of this program is used as input to the Histogram and Summary graph 
aragrams. 

■C Programmer : Keith J. Lirette 

Implementation of Musa - Basic Execution Time Model > 

=R0GRAM BASIC (INPUT, OUTPUT, in_file, out_file); 


CONST 

epsilon = 0.0*1; 

max_ints = 200; { max number of time intervals > 


TYPE 

array_type = ARRAYC1. .max_ints] of REAL; 
matrix_type = ARRAYC1..10, l..max_mts] of REAL; 

VAR 


Me, 

gamma_0, 

repeats, 

1 , 

k, 

runs, 
start , 

Me_of f set , 
num_pts, 

num_errors : INTEGER; 


{ number of bugs removed 3- 
{ number of total bugs > 

•C number of data generations per run > 

•C loop counters > 

{ number of runs desired > 

{ number of errors to 1st generate 3 
•C amt to increase # errors generated by > 
{ number of different # errors wanted 3- 


temp_sum, 

phi 

b0_est , 
bl_est , 
te, 

time_sum 

t 

in_f lie, 
out file 


{ dummy variable } 

REAL; i failure rate per bug, input J 

{ estimate of bO } 

■C estimate of bl } 

{ time to end of testing > 

matrix_type; 

array_type; { total time to failures > 

TEXT; 


{ 

FUNCTION get_b0_est (Me : INTEGER; 

te, 

bl_est : REAL) : REAL; 

C function to get estimate of bO given bl_est > 
begin 

get_b0_est := Me / ( 1 - exp(-bl_est * te) ); 
end; 
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BAE .P 


r\ 


{ > 

FUNCTION calc_bl_est (Me : INTEGER; { # errors -found } 

te, { time to end of testing > 

time_sum : REAL) : REAL; { sum of failure times > 

C function to get estimate of bl given failure data 3- 

VAR 

bl_est , 
low, 
high , 
lhs_est , 
rhs_est , 

Me x te 


FUNCTION lhs (Me : INTEGER) : REAL; 

•C function to get value of left hand side of equation solving 
for bl > 


{ holds working estimate of bl } 

•C low end of estimation interval } 

C high end of estimation interval > 

{ var to hold lhs of eq value > 

{ var to hold rhs of eq value > 

: REAL; { var to hold product calculated once } 


begin 

lhs := Me; 
end; i lhs } 

{ y 

FUNCTION rhs (Me x te, 
te, 

time_sum, 

bl_est : REAL) : REAL; 

t function to get value of right hand side of equation solving 
for bl } 


begin 

rhs := bl_est * ( Me_x_te / (exp(bl_est * te) - 1) + time_sum); 
end; 

C > 

begin 

{ initial guess at bl to start calculations ) 
bl_est := 0.0001; 

■C calculate product once and store for use > 

Me_x_te := Me * te; 

■C get left hand side, right hand side guess values > 
lhs_est ;= lhs (Me); 

rhs_est := rhs (Me_x_te, te, time_sum, bl_est); 

•C if lhs_est > rhs_est increase bl_est until not true > 

while lhs_est > rhs_est do 

begin 

bl_est := bl_est + 0.0001; 
lhs_est ;= lhs (Me); 

rhs_est := rhs (Me_x_te, te, time_sum, bl_est); 
end; 
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p 


law := bl_est - 0.0001; 
high := bl_est; 

{ perform bisection until value at which relationship -first 
changes is found 

> 

repeat 

bl_est := (low + high) / 2; 
lhs_est := lhs (Me); 

rhs_est := rhs (Me_x_te, te, time_sum, bl_est); 
if lhs_est < rhs_est then 
high := bl_est 
else 

low := bl_est; 

until abs(lhs_est - rhs_est) <= epsilon; 

calc_bl_est := bl_est; 
end; 


r 

\ 

i procedure to generate a set of failure times using number of repeats 
indicated. 


s 


procedure new_data ( gamma_0, 

Me, 
repeats 
var time_sum 
phi 

var t 


INTEGER; 

REAL; 

REAL; 

array_type) ; 


var 

i, j : INTEGER; 

previous : REAL; 

t_avg : array_type; 

bust : boolean; 


{ indices > 

{ value of previous time > 

{ array to hold running total to be averaged 
{ indicates overflow condition > 


begin 

•C initialize array that will hold running totals of times > 
for j := 1 to Me do 
t_avg[j3 := 0.0; 

J := 0; 

I get number of time samples indicated by variable repeats 3 
repeat 

{ initialize variables } 
bust := false; 
previous := 0.0; 
i := 0; 

•C generate required times > 
repeat 

i := l + 1; 

•C test for overflow condition > 
if (previous*phi > 719.0) then 
begin 

bust := true; 

end 
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else { generate next time interval > 
begin 

tCi] := previous + (1.0/phi) * 

In ( 1 - In (1-random (1.0) ) / (gamma_0 * exp (-phi 
previous := tCil; 
end; 

until (1 = Me) or bust; -C until all times gotten or over-flow } 

{. add time intervals to t_avg to keep running totals > 

if not bust then 

begin 

for 1 := 1 to Me do 

t_avgCi] := t_avg£i] + tCil; 

j := j + 1; 

end; 

until j = repeats; { repeat until required samples performed } 

£ calculate average times and total of time averages } 
time_sum := 0.0; 
for j := 1 to Me do 
begin 

tCj] := t_avg[j] / repeats; 
time_sum := time_sum + tCj]; 
end; 
nd; 


previous) ) ) 
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BEGIN 

{ prepare input and output file 3- 
reset (in_file) ; 
rewrite (out_f 1 le) ; 

while not eof(in_file) do 
begin 

C get parameters to produce output for > 
readln (in_f lie, gamma_0, phi, Me, repeats, runs); 

-C start point is 107. of total number of errors wanted > 
start := round ( (gamma_0 * 0.1) ); 

•C generate 107. more errors until total number of errors generated > 
Me_offset := round ( (gamma_0 * 0.1) ); 

{ calculate total number of error sets } 
num_pts := (Me - Start) div Me_offset + 1; 

C produce fault data for number of runs wanted } 

for j := 1 to runs do 

begin 

new_data (gamma_0, Me, repeats, temp_sum, phi, t); 

{ use subset of times generated to get estimates using 10, 20... 7. } 

for k := 1 to num_pts do 

begin 

{ calculate number of error times to use > 
num_errors := start + (k-1) * Me_offset; 

•C calculate time sum for errors used > 

time_sumCk, j3 := 0.0; 

for 1 := 1 to num_errors do 

time_sumCk, j3 := time_sumCk , j ] + tC13; 

•C set value for total testing time > 
teCk,j] := tCnum_errors3; 

•C get estimates for bl and bO > 

bl_estCk,j3 := calc_bl_est (num_errors, telk,j3, time_sumCk, j3) ; 
b0_estCk,j] := get_b0_est (num_errors, te[k,j], bl_esttk, j]) ; 
end; 

end; 

■C write outs set of calculated values > 

for k := 1 to num_pts do 

begin 

num_errors := start + (k-1) * Me_offset; 

writeln (out_f lie, gamma_0, ' ',phi,' num_errors,' ', repeats,' 
for j := 1 to runs do 

writeln (out_f lie, bl_estCk,j], b0_estCk,j3, time_sumCk, j3 , 


' , runs) 


end; 


teCk, j3) ; 







CONr .DENCE.P 


This program produces Confidence Interval summary graphs. The program reads input from a text 
file. Inputs are the total program faults (gamma_0) , failure intensity per fault (phi), the 
total number of faults to generate (He), and the number of runs to perform for each failure 
set size (runs). The program produces the desired number of runs starting with 10V. of the 
total program faults, increases by 107. until Me fault trials are runs. Confidence intervals 
(50, 55, ...,957.) are calculated for each of the trial sizes generated and summary graphs are 

plotted. 

•C Programmer : Keith 3. Lirette 

PROGRAM CONFIDENCE (INPUT, OUTPUT, in_gph, out_gph); 

CONST 

epsilon = 0.01; 
epsilonC = 0.00000001; 
max_ints = 200; 
num_ci = 10; 
line_length = SO; 
markl = 

mark2 = 

combojnark = 
num_grps = 20; 

TYPE 

array _type = ARRAY! 1 .. max_ints] of REAL; 
matrix_type = ARRAY! 1.. 10, l..max_ints3 of REAL; 
line_type = ARRAY!1. . Iine_length3 of char; 

k_rec_type = record 

ci : integer; 
val : real; 
end; 

k_arr_type = ARRAY! 1 . . num_ci 3 of k_rec_type; 

grp_rec_type = record 

value : real; 
labell ; integer; 
end; 

data_rec_type = record 

b0_in_pct, 
bl_in_pct , 
bl_low_mean, 
bl_high_mean, 
b0_mean , 
bl_mean, 
bl_high, 
bi_low, 
b0_min_mean, 
b0_max_mean ; real; 
labell : line_type; 

used_l , 

used_2 : boolean; 

end; 

grp_array_type = array! 1 .. num_grps3 of grp_rec_type; 
data_array_type = array! 1 .. num_grps3 of data_rec_type; 


C max number of time intervals } 

{ number of confidence intervals } 


{ max number of percentage divisions } 
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VAR 

Me, 

gamma_0 
time_sum, 
phi , 

I_val , 
bl_of f set , 
bO_min , 
bO_max , 
bl_min , 
bl_max , 
bl_mean , 
bl_min_sum, 
bl_max_sum, 
bO_min_sum, 
bO_max_sum, 
bl_sum, 
bO_sum 

repeats, 

ii. 

J > 

bl_in_cnt , 
bO_in_cnt , 
bO_untest, 
bO_min_lower , 

bl_high , 
bl_low, 
runs, 
start, 

Me_of f set , 
num_pts, 

1 , 

num_errors, 

k 

k val 


bO_est , 
bl_est, 
te 

pos_array 


data 

t 

line 

in_gph, 

out_gph 




CON 3ENCE.P 


£ number of bugs removed 3 
: INTEGER; £ number o-f total bugs 3 

{ sum o-f all -failure times } 

£ -failure rate per bug, input > 

£ -fisher in-formation o-f bl_est } 

£ magnitude of confidence interval around bl } 
£ lower value of confidence interval 3 
£ upper value of confidence interval 3 
£ lower value of confidence interval 3 
£ upper value of confidence interval 3 
£ average of bl over number of runs 3 
£ sum variable used to calc bl min mean > 

£ sum variable used to calc bl max mean 3 

£ sum variable used to calc bO min mean } 

£ sum variable used to calc bO max mean > 

£ sum variable used to calc bl mean 3 

: REAL; £ sum variable used to calc bO mean } 

C number of data generations per run > 


£ # of bl ests in own conf interval > 

£ # of bO ests in own conf interval > 

£ # of bO with no upper bound 3 
£ # of bO with no upper bound but 

whose lower bound is < estimate 3 
£ # bl ests above own conf interval > 

£ # bl ests below own conf interval 3 

£ number of runs desired 3 
£ number of errors to generate first 3 
£ amt to increase number of errors by 3 
£ number of different # of errors 1 
C counter } 

£ nummber of errors generated } 

: INTEGER; £ counter > 

: k_arr_type; £ holds confidence intervals and 

corresponding k values 3- 


£ estimate of bO 3 
£ estimate of bl 3 

: matrix_type; £ time to end of testing 3 

: grp_array_type; £ holds % divisions for graph and 

labels for y axis 3 


data_array_type; 
array _type; 
lme_type; 

TEXT; 


£ holds calculated fault data and 
graph labels 3 
£ total time to failures 3 

£ graph line built to print 3 
£ text file of input parameters 3 
£ output file for generated graphs 3 
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{ 

FUNCTION I (Me : INTEGER; 

bl_est, 

te” : REAL) : REAL; 


{ -function to get fisher information for bl } 


begin 

I := Me * ( 
end; 


1 / SQR (bl_est) - 

SQR(te) * exp(bl_est * te) / SQR( exp(bl_est * te) - 1) ); 


•C } 

FUNCTION get_bO_est (Me : INTEGER; 

te,. 

bl_est : REAL) : REAL; 

{ function to get estimate of bO given bl_est > 


begin 

get_bO_est := Me / ( 1 - exp (-bl_est * te) ); 
end; 


FUNCTION calc_bl_est (Me : INTEGER; 

te, 

time_sum : REAL) : REAL; 


{ # errors found } 
i time to end of testing } 
{ sum of failure times > 


{ function to get estimate of bl given failure data > 


VAR 

bl_est, 
low, 
high , 
lhs_est, 
rhs_est , 
Me x te 


C holds working estimate of bl } 

{ low end of estimation interval > 

{ high end of estimation interval } 

{ var to hold lhs of eq value > 

{ var to hold rhs of eq value } 

REAL; { var to hold product calculated once 


> 


{ } 

FUNCTION lhs (Me ; INTEGER) ; REAL; 


•C function to get value of left hand side of equation solving 
for bl }• 

begin 

lhs := Me; 
end; i lhs } 

* > 

FUNCTION rhs(Me_x_te, 
te, 

time_sum, 

bl.est : REAL) ; REAL; 


C function to get value of right hand side of equation solving 
for bl } 

begin 

rhs := bl_est * ( Me_x_te / (exp(bl_est * te) - 1) + time_sum); 
end; 
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{ } 

begin 

£ initial guess at bl to start calculations } 
bl_est := 0.0001; 

£ calculate product once and store -for use } 

Me_x_te := Me * te; 

£ get le-ft hand side, right hand side guess values } 
lhs_est := lhs (Me); 

rhs_e5t := rhs (Me_x_te, te, time_sum, bl_est); 

£ if lhs_est > rhs_est increase bl_est until not true > 

while lhs_est > rhs_est do 

begin 

bl_est := bl_est + 0.0001; 
lhs_est := lhs (Me); 

rh5_est := rhs (Me_x_te, te, time_sum, bl_est); 

end; 

low := bl_est - 0.0001; 
high := bl_est; 

£ perform bisection until value at which relationship first 
changes is found 

} 

repeat 

bl_est := (low + high) / 2; 
lhs_est := lhs (Me); 

rh5_est ;= rhs (Me_x_te, te, time_sum, bl_est); 

if lhs_est < rhs_est then 
high := bl_est 
else 

low ;= bl_est; 

until abs(lhs_est - rhs_est) <= epsilon; 
calc_bl_est := bl_est; 
end; 
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procedure new_data 


var 


gamma_0, 

Me : INTEGER; 
var time_sum : REAL; 

phi : REAL; 

var t : array_type); 


1 : INTEGER; { index } 

previous : REAL; { value of previous time generated } 

busted : boolean; { indicates if overflow condition present > 

begin 

{ repeat until valid times generated 3 
repeat 

•C initialize values } 
previous := 0.0; 
time_sum := 0.0; 


busted := false; 
l : = 0; 


•C generate times until done or overflow condition > 

while (i < Me) and not busted do 

begin 


i 


+ 1 ; 


C calculate new time T 

till := previous + (1.0/phi) * 

In ( 1 - In (1-random (1.0) ) / (gamma_0 * exp (-phi * previous) ) ) 
previous := tCil; 

{ add new time value to running total } 
time_sum := time_sum + tCi3; 


{ test for overflow condition } 
busted := ( phi * previous) > 88.0; 
end; 

until not busted; 
end; 


{ set confidence interval sizes and corresponing k values 3 
procedure initial ize_k_val (var k_val : k_arr_type) ; 
begin 

k_valC13.ci := 50; 
k_val C 1 3 . val ;= 0.675; 
k_val [23. ci := 55; 
k_valC21.val := 0.755; 
k_val[33.ci ;= 60; 
k_valC33.val ;= 0.84; 
k_valC43.ci ;= 65; 
k_valC41.val := 0.935; 
k_valC53.ci ;= 70; 
k_valC53.val ;= 1.035; 
k_valC63.ci := 75; 
k_valC6].val ;= 1.15; 
k valC73.ci := 80; 


} 
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k_val C73.val := 1.28; 
k_valC83.ci := 85; 
k_valC8].val ;= 1.44; 
k_valC93.ci := 90; 
k_valC93.val ;= 1.645; 
k_valC103.ci ;= 95; 
k_valC103.val ;= 1.96; 
and; 

{ > 

{ initialize array to hold 7. group break out and y-axis labels } 
procedure imtial_pos_array (var pos_array ! grp_array_type) ; 

var l ; integer; 

begin 

po5_arrayC13. value ;= 0.999; 
pos_arrayC13. label 1 := 100; 
for l ;= 2 to num_grps do 
begin 

pos_arrayCi3.value := pas_arrayCi-l 3 . value - 0.05; 
pos_arrayti 3. labell := pos_arrayCi-l 3 . labell - 5; 
end ; 
end; 

r \ 

t ' “ j 

{ procedure to convert integer < 1000 to a string } 
procedure strx (tag : integer; var val_str : line_type); 

var temp : integer; 

begin 

temp ;= tag div 100; 
val_strC13 ;= chr (48+temp) ; 
tag := tag - temp*100; 
temp := tag div 10; 
val_strE23 ;= chr (48+temp) ; 
tag := tag - temp*10; 
temp ;= tag; 

val_strC33 := chr (48+temp) ; 
end; 
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{ > 

{ procedure to get data line to print on graph } 
procedure get_line( pos_array : grp_array_type; 

var data : data_array_type; 
grp_num, 

data_pts : integer; 
highlight : boolean); 

var 

i , 

mark_pos : INTEGER; 
line : arrayE1..503 of char; 
fill : CHAR; 
begin 

•C highlight indicates this line is the line matching the confidence 
. interval size, so mark it with dashes > 
if highlight then 
fill := 
else 

fill := ' '; 

C initialize line > 
for i := 1 to data_pts*5 do 
1 meCi 3 := f ill ; 

•C if line not on confidence interval size, set dots every 5 so 
reading is easier > 
if not highlight then 

for l := 1 to data_pts do 
1 ineti*5 - 23 := ' . ' ; 

write <out_gph, pos_arrayCgrp_num3.1abell:4, ' !'); 

C determine if any data points fall in indicated X group and mark 3- 

far l := 1 to data_pts do 

begin 

mark_pos := i*5 - 2; 

if (pos_arrayCgrp_num3. value < dataCi 3.b0_in_pct) then 

if not dataCi 3. used_l then 
begin 

1 ineCmark_pos3 := markl; 
dataCi 3. used_l := true; 
end; 

if (pos_arrayCgrp_num3. value < datati3.bl_in_pct) then 

if not dataCi 3. used_2 then 
begin 

if lineCmark_pos3 = markl then 
lineCmark_pos3 := combo_mark 
else 

line[mark_pos3 := mark2; 
dataCi 3 . used_2 := true; 
end; 


end; 
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{ print line to output file > 
for 1 := 1 to data_pts*5 do 

write (out_gph, lineCil); 
writeln (out_gph) ; 
end; 

C > 

BEGIN 

reset (in_gph) ; 
rewr 1 te (out_gph ) ; 

imtial_pos_array (pos_array) ; 
initialize_k_val (k_val ) ; 

while not eof(in_gph) do 
begin 

readln (in_gph, gamma_0, phi, Me, runs); 

{ determine number of errors to first generate > 
start ;= round! (gamma_0 * 0.1) ); 

{ determine amt to increase errors by } 

Me_offset := round! (gamma_0 * 0.1) ); 

•C determine number of different errors set sizes > 
num_pts := (Me - Start) div Me_offset + 1; 

•C generate data for number of desired runs — total errors } 

for j ;= 1 to runs do 

begin 

new_data (gamma_0, Me, time_sum, phi, t); 

{ for each error set size, use # of errors from those errors 
generated and calculate bO and bl estimates > 
for k := 1 to num_pts do 
begin 

num_errors := start + (k-1) * Me_offset; 
time_sum := 0.0; 
for 1 := 1 to num_errors do 
time_sum ;= time_sum + till; 
teCk, j] ;= tCnum_errorsl; 

bl_estCk, jl ;= calc_bl_est !num_errors, teCk, j 1 , time_sum); 
b0_esttk, j] ;= get_b0_est (num_errors , teCk, jl, bl_estCk, jl); 
end; 
end; 

{ calculate averages for generated estimates } 

for k := 1 to num_pts do 

begin 

bl_sum := 0.0; 

b0_sum ;= 0.0; 

for j := 1 to runs da 

begin 

bl_sum ;= bl_sum + bl_est[k,j3; 
b0_sum := b0_sum + b0_estCk,j3; 
end; 

dataCk3.bl_mean := bl_sum / runs; 
dataCkl.b0_mean := b0_sum / runs; 
end; 
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{ test data on different cc intervals in k_val } 

for ii := 1 to num_ci do 

begin 

{ for each error set size, generate conf intervals on values } 

for k := 1 to num_pts do 

begin 

num_errors := start + (k— 1 ) * Me_offset; 

bO_in_cnt : = 0; 

bO_untest := 0; 

bO_min_lower := 0; 

bl_in_cnt := 0; 

bl_high := 0; 

bl_low : - 0; 

bl_min_sum : = 0.0; 

bl_max_sum : = 0.0; 

b0_min_sum : = 0.0; 

b0_max_sum := 0.0; 

{ for each data run, calculate confidence and test } 

for j := 1 to runs do 

begin 

I_val := I (num_errors, bl_est[k,j], te[k,j]); 
bl_offset := k_val[ii] . val / SQRT(I_val); 
bl_min := bl_est[k, j] - bl_offset; 
bl_max := bl_est[k, j] + bl_offset; 
bl_min_sum := bl_min_sum + bl_min; 
bl_max_sum bl_max_sum + bl_max; 
if bl_min <= epsilon2 then 
begin 

bl_min := 0.0; 

b0_untest := b0_untest + 1; 

if get_b0_est(num_errors, te[k, j], bl_max) < gamma_0 then 
b0_min_lower := b0_min_lower + 1; 

end 

else 

begin 

b0_min := get_b0_est(num_errors, te[k, j], bl_max); 
b0_min_sum := b0_min_sum + b0_min; 
b0_max := get_b0_est(num_errors, te[k, j], bl_min); 
b0_max_sum b0_max_sum + b0_max; 

if (gamma_0 >= b0_min) and 

(gamma_0 <= bO_max) then 
b0_in_cnt := b0_in_cnt + 1; 

end; 

if phi < bl_min then 

bl_high := bl_high + 1 
else if phi > bl_max then 
bl_low := bl_low + 1 
else 

bl_in_cnt := bl_in_cnt + 1; 


end; 

strx(num_errors, data[k] .labell); 

data[k] .b0_in_pct := (b0_in_cnt + b0_min_lower)/runs; 
data[k] .bl_in_pct := bl_in_cnt / runs; 
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dataCk3.used_l := false; 
dataCk]. used_2 := false; 
dataCk]. bl_high := bl_high / runs; 
dataCk] . bl_low := bl_low / runs; 

dataCk]. bl_low_mean := bl_min_sum / runs; 
dataCk]. bl_high_mean := bl_max_sum / runs; 

if bO_untest < runs then 
begin 

dataCk] .bO_min_mean := bO_min_sum / (runs - bO_untest); 
if bO_untest >= (runs div 2) then 
dataCk ].bO_max_mean := 9999.0 
else 

dataCk] . bO_max_mean := bO_max_sum / (runs - bO_untest); 

end 

else 

begin 

dataCk]. bO_min_mean := 0.0; 
dataCk]. b0_max_mean := 9999.0; 
end ; 

end; 

C print results for one conf interval size > 

writeln (out_gph, ' Graph of 7. of bO/bl Estimates Falling in 

k_val Ci l ] . ci : 3, '7. Cl '); 

writeln (out_gph) ; 

writeln(out_gph, ’ # => bO @ => bl * => bO ?< b 1 ' ) ; 

writeln (out_gph) ; 

writeln (out_gph, ' Gamma_0 = ' ,gamma_0:3, ' Phi = *,phi:5:4, ' ':5, 

'runs : ',runs:2); 
writeln (out_gph) ; 

writeln (out_gph) ; 

writeln (out_gph, ' 7. in'); 

writeln (out_gph) ; 

for j := 1 to num_grps do 

begin 

get_l ine(pos_array, data, j, num_pts, k_valCn].ci = pos_arrayCj].labell) ; 
end; 

write (out_gph , ' 

for k := 1 to num_pts do 


write (out_gph , ' ') ; 

writeln (out_gph) ; 

for 1 := 1 to 3 do 

begin i 


write (out_gph, ' ':6); 
for k := 1 to num_pts do 

write (out_gph , ' dataCk]. label 1 Cl ], ' '); 

if 1 = 2 then 

writeln(out_gph, ' # faults fixed') 

else 

writeln (out_gph) ; 

end; 

writeln (out_gph) ; 

writeln (out_gph , ' Mean 7.bl 


G-15 



CONr *i)ENCE.P 


*/.bO 

writeln(out_gph, ' Me bl 

' low in ci 
writeln (out_gph , ' 


Mean ' 

bl Min bl Max 
bO bO min 


); 

in ci 
bO max ' ) 


') 


■for k := 1 to num_pts do 

writeln (out_gph , dataCk] . labell , ' dataCk]. bl_mean:7:6, ' 
dataCk]. bl_low_mean:7:6, ' ' , dataCk]. bl_high_mean:7:£j, ' 

dataCk]. bl_m_pct: 4: 3, ' 

dataCk]. bl_high:4:3, ' dataCk]. bl_low:4:3, ' 
dataCk]. b0_in_pct:4:3, ' ' , dataCk]. b0_mean:6: 1, ' 

dataCk].bO_min_mean:6: 1, ' ' ,dataCk].bO_max_mean:8: 1) ; 

wri teln (out_gph ,chr ( 12) ) ; 
end; 
end; 

close (in_gph) ; 
close(out_gph) ; 

END. 


high ' , 

9 


1 
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This program performs no data generation. It produces summary graphs of the of data generated 
by the BASIC. P program. 

PROGRAM GRAPH; 

const 

line_length = 65; 
mark = ; 

num_grps = 31; 
max_pts = 20; 

type 

array_type = array! 1 . . 1003 of REAL; 

Iine_type = STRINGElme_length]; 

grp_rec_type = record 

value : real; 
labell : integer; 
end; 

data_rec_type = record 

value, 
bO_mean , 
bO_std , 
b l_mean , 
bl_std, 
tsum_mean , 
te_mean : real ; 
labell : stringC33; 
used : boolean; 
end; 

grp_array_type = arrayC 1 . . num_grps] of grp_rec_type; 

data_array_type = array [ 1 .. max_ptsl of data_rec_type; 


<■ number horizontal graph lines > 

{ max number of vertical graph lines } 


var 

gamma_0, 

repeats, 

runs, 

Me, 

i , 

^ » 

temp_length, 

j " : INTEGER; 

phi : REAL; 

pos_array : grp_array_type; 

data : data array type; 

bl, 

te, 

tsum, 

bO : array_type; 

line : line_type; 

in_f lie : TEXT? 

out_f lie : TEXT; 

f ile_name, 

file_name2 : STRINGC143; 

temp_str : STRINGC33; 
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{ > 

{ This array holds the break out percentages -for graph lines and the 
corresponding y-axis labels 

} 

procedure imtial_pos_array (var pos_array : grp_array_type) ; 
begin 

pos_arrayCl]. value := 1.95; 
pos_arrayC13.1abell : = 100; 

•for 1 := 2 to 10 do 
begin 

pos_array[i 3. value := pos_array£i-13. value - 0.1; 
pos_arrayCi 3. labell := pos_arrayCi-13. labell - 10; 
end; 

pos_arrayCll].value := 1.045; 
pos_arrayC113.1abell := 5; 

■for i := 12 to 21 do 
begin 

pos_array[i3.value := pos_arrayCi-13. value - 0.01; 
pos_arrayCi 3 . label 1 := pos_arrayCi-13.1abell - 1; 
end; 

pos_arrayC223. value := 0.85; 
pos_array[223 . label 1 := -10; 

•for l := 23 to 30 do 
begin 

pos_arrayCi 3. value := pos_arrayti-13. value - 0.1; 
po5_arrayCi 3. labell := pos_array[i-13.1abell - 10; 
end; 

pos_arrayC313.value := 0.0; 
pos_arrayC31 3 . labell := -100; 
end; 

C > 


FUNCTION Mean (Values: array_type; Count: integer) : Real; 

(* This FUNCTION calculates the MEAN of count number of VALUES *) 


Var 

Accum : Real; (* sum of all values *) 

I : Integer; (* loop control variable *) 

Begin (* mean *) 

Accum := 0.0; (* initialize sum *) 

for I := 1 to Count do 

Accum := Accum + ValuesCI3; 

Mean := Accum / Count; (* calculate MEAN *) 

End; (* mean *) 

■C > 

FUNCTION Standev (Values: array_type; Mean: REAL; Count: INTEGER) : REAL; 

(* This FUNCTION calulates the Standard Deviation of COUNT 
amount VALUES gives the MEAN *) 


Var 

Sqr_Accum, 
temp : real; 

I : Integer; 


(* sum of all values squared *) 
(* loop control variable *) 
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Begin (* Standev *) 

Sqr_Accum := 0.0; (* Reset to 0 #) 

For I := 1 To Count Do (* sum up squared values *) 

begin 

temp := (Values! 1 3 - Mean); 

Sqr_Accum := Sqr_Accum + temp * temp; 
end; 

Standev := SORT < Sqr_Accum / Count) (* Figure Standard Deviation *) 

End; (* Standev *) 

{ > 

{ This procedure prints a horizontal graph line after all data has been 
totalled 

J 

procedure get_line( pos_array : grp_array_type; 

var data : data_array_type; 
grp_num, 

data_pts : integer); 

var 

i > 

mark_pos : INTEGER; 
line : STRINGCS03; 

■fill : CHAR; 
begin 

if pos arrayCgrp numl.labell = 0 then 
■Fill ~ c = 
else 

fill := ' '; 

for l := 1 to data_pts*5 do 
lineCi 3 := fill; 

write (out _f l le,pos_arrayCgrp_num3. label 1:4,' ! ' ) ; 

for i := 1 to data_pts do 

begin 

if (pos_arrayCgrp_num3.value < dataCi 3. value) then 
if not dataCi 3. used then 
begin 

lineCi*5 - 13 := mark; 
dataCi 3. used := true; 
end 
else 

line[i*5 - 13 := 

end; 

for l := 1 to data_pts*5 do 
write(out_file, lineCi3); 
writeln(out_file) ; 
end; 

i } 

begin 

imtial_pos_array (pos_array) ; 

write ('Enter Name of Data File => '); 

read 1 n ( f i 1 e_name ) ; 

assign (in_f lie, file_name); 

reset (in_f lie) ; 

f i le_name2 : = ' ' ; 
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write( 'Enter Name -for Output File (Ret -for .gph ext) => '); 

readln (f 1 1 e_name2) ; 

if length (fi le_name2) = 0 then 

file_name2 := concat (f 1 le_name, ' . gph ' ) ; 
assign (out_fi le, file_name2); 
rewrite (out_f lie) ; 
k := 0; 

while not eof(in_file) do 
begin 

k := k + 1; 

readln (in_fi le, gamma_0, phi, Me, repeats, runs); 

for I := 1 to runs do 

begin 

readln (in_f lie, blEil, bOCi], tsumCi], teCi]); 
end; 

str (repeats, temp_str); 
temp_length := length (temp_str) ; 

dataCk] . label 1 := concat (copy ( ' ',1,3 - temp_length) , temp_str) ; 

dataCk] .b0_mean := MEAN(bO, runs); 

dataCk]. bljnean := MEAN(bl, runs); 

dataCk] . tsum_mean := MEAN(tsum, runs); 

dataCk ].te_mean := MEAN(te, runs); 

dataCk]. b0_std := Standev(bO, dataCk] . bO_mean , runs); 
dataCk] . bl_std := Standev(bl, dataCk]. bl_mean, runs); 
dataCk] . value := dataCk] . b0_mean/gamma_0; 
dataCk]. used := false; 
end; 

writeln (aut_f lie, ' Graph of '/. bO estimate differs from gammaO : ') 
writeln (out_f lie) ; 

writeln (out_f lie, ' Gamma_0 = ',gamma_0:3, ' Phi = ',phi:5:3,' ';10, 
'Source : ' ,f l le_name) ; 
writeln (out_f lie, ' n = ',Me:3,' ':23, 

'Output ; ' ,f ile_name2) ; 

writeln (out_f lie) ; 

writeln (out_fi le, ' '/.off'); 

writeln (out_f lie) ; 

for i := 1 to num_grps do 

begin 

get_line(pos_array, data, i, k); 
end; 

write(out_f lie, ' ':5,'-')j 
far l := 1 to k do 

write (out_f lie, ' ' ) ; 

writeln (out_f lie) ; 

for i ;= 1 to 3 do 
begin 

write(out_f lie, ' ':7); 
for j := 1 to k do 

write (out_fi le, ' ' ,dataCj]. labell Cl ] , ' '); 

if i = 2 then 

writeln (out_f lie, ' # repeats') 

else 

writeln (out_f i le) ; 

end; 
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writeln (out_fi le) ; 


{ print statistical results } 
■for 1 := 1 to 94 do 
write (out_f lie, ) ; 
wri teln (out_f l le) ; 
writeln (out_f lie, ' 

writeln (out_f lie, ' Me 

Std 

writeln (out_f lie, ' 


bl 

tsum te'); 

Mean Std 

Mean Mean ' ) ; 




bO' 


Mean 




for l := 1 to k do 

writeln (out_f lie, ' ':5,dataCi].labell , ' ', dataCi3.bl_mean:9:8, ' 

dataCi 3 . bl_std: 9:8, ’ ', dataCi3.bO_mean:7;2, ' 

dataCi 3. b0_std:7:2, ' dataCi 3. tsum_mean: 10: 1 , ' 

dataCi 3.te_mean:8: 1) ; 


close(in_f lie) ; 
close (out_f l 1 e) ; 
and. 
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This program generates scatter plots (histograms of the data generated by the BASIC. P program. 

PROGRAM HISTOGRAM; 
const 

mark = '*'; 
num_grps = 20; 

type 

table_type = array! 1 .. num_grps, 1 .. 100] of Char; 

count_array_type = array C 1. .num_grps] of INTEGER; 

/ar 

gamma_0, 
repeats , 
runs, 

Me, 

group , 

ma>:_grp_cnt , 
i » 

j : INTEGER; 


phi , 


bl_est , 


bl_max , 


bl_min , 


b0_est : 

REAL; 

data : 

TEXT; 

out_file : 

TEXT; 

f ile_name, 


file_name2 : 

STRING! 14] ; 

count_arr : 

count_array_type 

table : 

table type; 

graph_type : 

CHAR; 

b0_graph : 

r , 

BOOLEAN; 

procedure initialize; 

var 


l, j : INTEGER; 

begin 



for i := 1 to num_grps do 
begin 

count_arr!i] := 0; 
for j := 1 to 100 do 
tableCi,j3 := * '; 

end; 

end; 


begin 

write ('Enter Name of Data File => '); 
readln (f ile_name) ; 
assign (data, file_name); 
reset (data) ; 

write( 'Enter Name for Output File => '); 

readln (fi le_name2) ; 

write ( 'Graph for bO (Y/N) '); 
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readln (graph_type> ; 

bO_graph := graph_type in C'Y','y'3; 

assign (out_fi le, file_name2); 

rewrite(out_f lie) ; 

while not eof (data) do 

begin 

initialize; 

readln (data, gamma_0, phi, Me, repeats, runs); 

•for I := 1 to runs do 
begin 

readln (data, bl_est, bO_est); 
l-f bO_graph then 

group := 1 + trunc( ( (b0_est/gamma_0) * 100.0 - 5.0) / 10) 
else 

group := 1 + trunc( ( (bl_est/phi) * 100.0 - 5.0) / 10); 
if group > num_grps then 
group := num_grps 
else if group < 1 then 
group : = 1 ; 

count_arr [group! := count_arr [group] + 1; 
tableCgroup,count_arr[group]3 := mark; 
end; 

max_grp_cnt := 0; 

for I := 1 to num_grps do 

if count_arrCI] > max_grp_cnt then 
max_grp_cnt := count_arrCI]; 
if max_grp_cnt < 48 then 

max_grp_cnt := max_grp_cnt + 2; 
if b0_graph then 

writeln (out_f l le, ' Chart of Estimates of bO with : ’) 

else 

writeln (out_fi le, ' Chart of Estimates of bl with : '); 

writeln (out_f lie) ; 

writeln (out_f lie, ' Gamma_0 = ',gamma_0:3, ' Phi = ',phi:5:3,' ':10, 

'Source : ' ,f l le_name) ; 

writeln (out_f lie, ' n = ',Me:3, ' r = ', repeats: 3, ' ' : 12 , 

'Output File : ' ,f ile_name2) ; 
writeln (out_f lie) ; 
for J := 1 to max_grp_cnt do 
begin 

write(out_f lie, max_grp_cnt - J + 1:2,' .'); 
for I := 1 to num_grps do 

write(out_file, ' ', tablet I ,max_grp_cnt - J + 1],' '); 
writeln (out_f lie, ' .'); 
end; 

write(out_f lie, ' -'); 

for I := 1 to num_grps do 

wnte(out_f lie, ' ') ; 

writeln (out_f lie, ' — ') ; 

writeln (out’file, ' -90 80 70 60 50 40 30 20 10 0 10 20 30 40 50 60 70 80 90 100+') 

writeln (out_f l le) ; 

writeln (out_f lie, ' 7. Estimate Differs from Actual Value '); 

writeln (out_f ile,chr (12) ) ; 
end; 

close(data) ; 
close(out_file) ; 
end. 
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- Code/Math Cross-Ref 



Equations used in BASIC. P to generate Failure Data 


This function is an implementation of equation 7. 


Po 


1 - expCf^t,,) 


Here, 


Me - equals the number of faults found, 
te - equals the time to end of testing, 
bl_est - equals the calculated estimate for bl 


FUNCTION get_bO_est (Me : INTEGER; 

• te, 

bl_est : REAL) : REAL; 

{ function to get estimate of bO given bl_est ) 
begin 

get_bO_est := Me / ( 1 - exp (-bl_est * te) ); 
end; 


This function is an implementation of equation 11. 


Ell 

Pi 


* 3*0 

exp -1 



= 0 


The equation was rearranged to the following form: 


= Pi 


exp (Pit,,) -1 M 




The function has 2 subfunctions, lhs and rhs, which represent each side of the 
rearranged equation. The bisection method is used to find a point where the 
values of the lhs and rhs functions are less than a given epsilon. 

FUNCTION calc_bl_est (Me : INTEGER; { # errors found } 

te, { time to end of testing ) 

time_sum : REAL) : REAL; { sum of failure times } 

{ function to get estimate of bl given failure data } 

VAR 

bl_est, { holds working estimate of bl } 

low, { low end of estimation interval } 

high, { high end of estimation interval } 

lhs_est, { var to hold lhs of eq value } 

rhs_est, { var to hold rhs of eq value } 

Me_x_te : REAL; { var to hold product calculated once } 
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t } 

FUNCTION lhs (Me : INTEGER) : REAL; 


{ function to get value of left hand side of equation solving 
for bl } 

begin 

lhs : = Me ; 
end; { lhs } 


1 ) 

FUNCTION rhs(Me_X_te, 
te, 

time_sum, 

bl est : REAL) : REAL; 


{ function to get value of right hand side of equation solving 
for bl ) 


begin 

rhs := bl_est * ( Me_x_te / (exp(bl_est * te) - 1) + time_sum) ; 

end; 

( ) 

begin 

{ initial guess at bl to start calculations ) 
bl_est := 0.0001; 

{ calculate product once and store for use } 

Me_x_te := Me * te; 

{ get left hand side, right hand side guess values } 
lhs_est := lhs (Me) ; 

rhs_est := rhs (Me_x_te, te, time_sum, bl_est) ; 

( if lhs_est > rhs_est increase bl_est until not true } 

while lhs_est > rhs_est do 

begin 

bl_est := bl_est + 0.0001; 
lhs_est := lhs (Me) ; 

rhs_est := rhs (Me_x_te, te, time_sum, bl_est) ; 
end; 

low := bl_est - 0.0001; 
high := bl_est; 

{ perform bisection until value at which relationship first 
changes is found 

} 

repeat 

bl_est := (low + high) / 2; 
lhs_est := lhs (Me); 

rhs_est rhs (Me_x_te, te, time_sum, bl_est) ; 
if lhs_est < rhs_est then 
high : = bl_est 
else 

low := bl_est; 

until abs(lhs_est - rhs_est) <= epsilon; 

calc_bl_est := bl_est; 
end; 
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This procedure is an implentation of equation 23. The produces 'repeats' number 
of failure time data sets and averages the corresponding failure times. 


t = 



+ 


ln(l - z) 

v o-do/To ) tsum 

Yo e 


procedure new_data ( gamma_0. 

He, 
repeats 
var time_sum 
phi 

var t 


{ total number of faults } 

{ number faults to generate) 
INTEGER; { number of repetitions } 
REAL; { sum of all failure times } 

REAL; { failure rate per fault } 

array_type) ; { array of generated times } 


{ procedure to generate a set of failure times using number of repeats 

indicated. 

) 


var 

i, j : INTEGER; 

tsum : REAL; 

t_avg : array_type; 

bust : boolean; 


{ indices } 

{ value of previous failure time } 

{ array to hold running total to be averaged } 
{ indicates overflow condition } 


begin 

{ initialize array that will hold running totals of times } 
for j := 1 to Me do 
t_avg [ j ] : = 0.0; 

3 := 0 ; 

{ get number of time samples indicated by variable repeats } 
repeat 

{ initialize variables } 
bust := false; 
tsum := 0.0; 

i := 0; 


{ generate required times } 
repeat 

i := i + 1; 


{ test for overflow condition } 
if (tsum * phi > 719.0) then 
begin 

bust := true; 

end 


else { generate next time interval } 
begin 

{ 1.0/ phi = gamma_0 / lambda_0 } 
t[i] := tsum + (1.0/phi) * 

In ( 1 - In (1-random (1. 0) ) / 

(gamma_0 * exp (-phi * tsum) ) ) ; 
tsum := t [i] ; 
end; 

until (i = Me) or bust; ( until all times gotten or overflow } 

{ add time intervals to t_avg to keep running totals } 

if not bust then 

begin 

for i := 1 to Me do 


H-3 



t_avg[i] := t_avg[i] + t [i] ; 
j := j + 1; 
end; 

until j = repeats; { repeat until required samples performed } 

{ calculate average times and total of time averages } 
time_sum := 0.0; 
for j := 1 to Me do 
begin 

t[j] := t_avg [ j ] / repeats; 
time_sum ;= time_sum + t [ j ] ; 
end; 
end; 
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Additional Confidence Interval Program Equations 


This function is an implementation of equation 15, 


r<Pi> 



tjexp 

[exp (P 1 t„) -1 ] 2 


FUNCTION I (Me : INTEGER; 

bl_est, 

te : REAL) : REAL; 

{ function to get fisher information for bl } 
begin 


I := Me * 
end; 


1 / SQR (bl_est) -• 

SQR(te) * exp(bl_est * te) / SQR( exp (bl_est * te) - 1) 


% 


The following code segment calculates confidence interval boundaries implementing 
equation 14. 



I_val := I (num_errors, bl_est[k,j], te [k, j]); 
bl_offset := k_val [ii] . val / SQRT(I_val); 
bl_min := bl_est[k, j] - bl_offset; 
bl_max := bl_est[k, j] + bl_offset; 

bO_min := get_bO_est (num_errors, te[k, j], bl_max) ; 
bO_max := get_bO_est (num_errors, te[k, j], bl_min) ; 
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End of Document 



